p-value calculation method Background
Trying to understand the way to calculate the p-value. Please exclude hypothesis testing here because in my understanding, hypothesis testing as one p-value use case is a different topic from calculating p-value itself. I believe mixing them making things complicated unnecessarily.
Suppose there is a Hard Disk Drive (HDD) type 0 called $H_0$. The MTTF (Mean Time to Failure) data has been collected and its mean $\mu=2.97$ and $\sigma=1.54$..
H0 = np.array([1] * 10000 + [2] * 25000 + [3] * 4000 + [4] * 3000 + [5] * 20000)
np.random.shuffle(H0)

u = np.mean(H0)
sd = np.std(H0)
print(f"mean {u:.2f} std {sd:.2f}")
---
mean 2.97 std 1.54


Then, run random sampling (sample size 100) from $H_0$.

*

*Take 100 samples and calculate the sample mean $\overline{\mu}$ of MTTF.

*Repeat the step 1 many many times.

*Plot the histogram of the $\overline{\mu}$
def sample_means(a, N = 10000):
    means = []
    for i in range(N):
        means.append(np.mean(np.random.choice(a=a, size=100, replace=False)))
    return np.array(means)

means = sample_means(H0, 100000)
plt.hist(means, bins=50)
plt.grid()

sm = np.mean(means)
ssd = np.std(means)
print(f"sampled: mean {sm:.2f} sd {ssd:.2f}")
---
sampled: mean 2.97 sd 0.15

Update
According to CLT, the SD of the ND of sampling means $\overline {\sigma} = \sigma / \sqrt {n} = 0.1545$ where n is sample size (100).

This approximates the Normal Distribution.
p = scipy.stats.norm(loc=sm, scale=ssd)

x = np.linspace(mu - 3*ssd, sm + 3*ssd, 1000)
plt.plot(x, p.pdf(x))
plt.grid()

print(p.cdf(np.inf))
---
1.0


Question
As far as I understood, p-value is to calculate the probability where the mean of the sample exceeds a certain value m (for right side one tailed).
$P(h \geqslant m \mid H_0)$
This probability is the area in ND where the mean is m or larger, because the ND approximates the multiple samplings from $H_0$.
Then the formula to get p-value is:
$P(h \geqslant m \mid H_0) = (1 - \mathtt{p.cdf(m)})$.
def p_value(m):
    return 1 - p.cdf(m)

Is this correct?
Update
Found a nice example with actual calculations.

*

*Hypothesis Testing in Finance: Concept and Example
When the sampling mean is $\overline {x}$, the Z value of the sample mean is below:

The p-value is the area in the ND where the $z > Z$.
So it looks my understanding seems correct. Although the formula is provided in the context of Hypothesis Testing, calculating p-value itself is achieved based on the data $H_0$, $\mu$ (mean of $H_0$), $\sigma $ (sd of $H_0$), and a sample mean $\overline { x }$. Apparently it stands without Hypothesis Testing.
I am still not sure why p-value calculation steps (solely how to get it) should bring in the discussion of what Hypothesis Testing is, just because p-value is used there. I only need to understand the steps to get p-value.

I think we should start with a standard definition of a p-value.  They will vary from book to book, but let us choose a simple one.
A p-value is the probability of obtaining results that are as extreme or more extreme than the observed results given that the null hypothesis is correct.
So a p-value is not a probability in the traditional sense.  If you would change your hypothesis, then the p-value changes.  It cannot exist in a manner that is separate from the null hypothesis.
You are mostly correct in that you are sampling data and finding an estimator.
For example, if you created 1,000,000 experiments of a sample size of 30 with a known variance of one, and a hypothesized population mean of zero, then you could find the distribution of the sampling means.  That is the sampling distribution of the sample means whose with 30 observations in each experiment.

However, what if we had calculated the medians instead?  We would get a different set of probabilities.  Our p-value is contingent on our stated hypothesis, here $\mu=0$; our model, which here is a simple normal model; and our choice of the estimator, or equivalently our choice of a loss function.  There are very legitimate reasons to use the median instead of the mean, though it is also uncommon for normal data.

But our concern is really our z-score.  After all, if we observe a sample mean of 12, its interpretation is very different if $\sigma^2=1$ and $\sigma^2=10000.$
So, the sampling distribution of our z-statistic conditions our sample mean so that it scales in a meaningful manner.  In an important sense, because we did already know the variance happens to be one, that was not necessary but that is not always the case.

Finally, if our null hypothesis is true, then the p-values should be uniformly distributed between zero and 100%.

That counter-intuitive result is true for any hypothesis test using p-values.  If you want to make a mistake no more than 1% of the time, then you want no more than 1% of the p-values to be significant in error.  Likewise, if someone else wants to make a mistake no more than 5% of the time, then they don't want more than 5% of the p-values to be significant in error.  That would go as well for people that want a 10% or 25% error rate.
The only way to make that work is for p-values to be uniformly distributed when the null is true.
Do note that the sampling method here is not vertical at the ends because it is a kernel density estimate and it misunderstands that there is no data supposed to be past zero or one.
Most simple problems already have an analytic solution.  However, if it does not, then the following steps will work.
Create the histogram implicit in your null hypothesis by repeated sampling.  It is important to know that the sample size is important.  If you change your sample size, then you should be changing your sampling distribution if we ignore a handful of degenerate cases.
If you are doing a two tailed test, then calculate region, $[\mu-A,\mu+A]$ where the probability of being in that region equals or just exceeds your cutoff probability.
Anything outside of that rejects your null.
If you are doing a left-tailed test, $\mu\le{k}$, then calculate the probability of being in the region to the left of $k$.  Do the same in the opposite direction for $\mu\ge{k}$ for the right tailed test.
rm(list = ls())


#library calls

library(ggplot2)
library(ggthemes)

observations<-30
experiments<-1000000

#Generate sample of 30x1,000,000 with center zero and variance of one
x<-matrix(rnorm(observations*experiments),nrow = observations)

#generate sample mean over sample X across experiments
sample_mean<-apply(x,2,mean)
sample_median<-apply(x,2,median)
rm(x)

#create z-score
z<-(sample_mean)*sqrt(observations)

#create p-values
p_value<-2*pnorm(-abs(z))

#create plot
plotting_data_frame<data.frame(z_score=z,p_value=p_value,sample_mean=sample_mean,sample_median=sample_median)

a<-ggplot(data = plotting_data_frame,aes(sample_mean))+theme_economist()+geom_density()+labs(title = "Simulation of Sampling Density of the Mean")
print(a)

a<-ggplot(data = plotting_data_frame)+theme_economist()
a<-a+geom_density(aes(sample_mean),color="Red")
a<-a+geom_density(aes(sample_median),color="Black")
a<-a+labs(x="Estimator",y="Density",title="Comparison of Sampling Distribution of Mean and Median",subtitle = "Sampling Distribution of the Mean (red), Sampling Distribution of the Median (Black)")
print(a)

a<-ggplot(data=plotting_data_frame,aes(z_score))+theme_economist()
a<-a+geom_density()
a<-a+labs(title = "Sampling Distribution of the Z-Score")
print(a)

a<-ggplot(data = plotting_data_frame,aes(p_value))+theme_economist()
a<-a+geom_density()

print(a)

The basic idea is that if you covered the entire sample space that is implied by your null hypothesis, then your choice of estimator should create the sampling distribution of that estimator, including the impact of any nuisance statistics like the variance.
Let us walk through the example of a solved problem, which is testing the location of a sample mean, when the variance is unknown for a single sample of size $n$ where the data is normally distributed.
If we use the sample mean as our estimator, which is reasonable, then our test statistic is usually denoted $t_n$ and uses Student's t distribution to calculate the probability of getting an observation as extreme or more extreme than what we observed.
Our formula is $$t_n=\frac{\bar{x}-\mu_0}{\frac{s}{\sqrt{n}}},$$ with $df=n-1$.  The first formula contains our estimate of the mean, $\bar{x}$, our hypothesized mean $\mu_0$, and our scaling $\frac{s}{\sqrt{n}}.$  Our second equation gives us our degrees of freedom, which is used in Student's distribution.
Student's t-distribution saves you from having to create millions of samples.
Also note that p-values tell you nothing if the null hypothesis is false.  A low p-value such as $p<.001$ only says that if the null is true, then less than one in one thousand samples will be as extreme or more extreme than what you observed.
In non-technical language, if we reject the null because the p-value is so low, that is the equivalent to saying "if the null is true, then our data is pretty wonky, so our null is probably not true."
If you reject the null hypothesis, then you are really saying that you trust the quality of your data and not the correctness of the stated hypothesis.
A simple way to think about a p-value, since the area in a probability distribution must sum to one, is that it is one minus the percentage of area taken up by the null hypothesis until it just touches the sample statistic.
If it is a two sided test, then you start at the hypothesized point and work your way outward, equally in both directions, until you just touch the observed value.
For a one sided test, you start at positive or negative infinity (or a finite bound in some cases) and keep filling up the area under the curve until you touch the observed point.  That is the area of the null.  The p-value is 1 minus the area of the null.
