How to prove statistically that mean jumps when a parameter crosses an unknown threshold value? We have a random numbers generator (it is actually a physical device standing outside of our building). By collecting of about 3000 numbers, generated by the device, we came to the conclusion that the mean of the generated numbers is positive (although it is close to zero). The statistical significance of this conclusion is quite good.
Recently we started to suspect that when temperatures outside are high, the device tends to generate slightly smaller numbers. Even more, we suspect that when temperature is above a threshold temperature, the mean of the generated numbers is negative (although we do not know what is exactly the value of this threshold temperature).
I also need to add, that we have a reason to believe, that the dependency of the mean of the generated numbers on the temperature is a step function. To make it more clear, we believe that when temperature is high enough two parts of the device disconnect from each other and the device switches to a new "regime" (associated with the generation of presumably smaller numbers).
So, my question is, how can we calculate the statistical significance of our hypothesis that there is a temperature triggered switch in the mean of the generated numbers?
We could, for example, find a value of the threshold temperature empirically such that the mean of the numbers, generated by the device, is negative and as low as possible and then we could calculate the statistical significance that the mean of these values is really negative (in the statistical sense). However, I guess that this approach is false, since we have manually chosen the split such that numbers corresponding to the high temperature are small, so it is not fair to prove if they have statistically significant negative mean.
In other words, would a threshold value of the temperature be given, we could just check if we have a statistical evidence to believe that values, corresponding to the high temperatures, have negative mean.
ADDED:
My question is not about physics. The "physical device" does not exists. I made the example just to make it less abstract and easier to grasp.

Just in case a Bayesian approach is of interest, there are some potential tools readily available in R or Python for use. The one probably most appropriate for your case is the mcp package--regression with multiple changepoints. Another popular one is bcp. But here I will choose a package Rbeast (https://github.com/zhaokg/Rbeast) most familiar to me (bcz I wrote it) to illustrate the ideas using the sample data from @user4422:
set.seed(1)
n_obs <- 4000
temperatures                <- rnorm(n_obs)
random_numbers              <- rnorm(n_obs)

temperature_threshold       <- 1.5
increase_if_above_threshold <- 0.5
random_numbers[temperatures>temperature_threshold]=random_numbers[temperatures>temperature_threshold]  + increase_if_above_threshold           

df <- data.frame(random_numbers, temperatures)
df <- df[order(df$temperature),]


library(Rbeast)
out = beast( 
            df$random_numbers, 
            season='none',        # only consider a trend component
            torder.minmax=c(0,0), # the piecewise trend is modelled as flat lines (i.e., the min and max polynomial orders are both 0)
            tcp.minmax  =c(0,1) # as in your case, the possible range of number of changepoints are 0 to 1.
           );

plot(out, vars=c("y","t","tcp"), main='Posterior probability of changepoint in the mean trend')

In the 'prob' subplot below is a curve of probability of changepoint occurrence. Apparently, peaks correspond to locations where changepoints occur.

Print out some addition posterior statistics: The posterior probability of observing a mean shift at 3748 is 0.845.
print(out)

#####################################################################
#                      Trend  Changepoints                          #
#####################################################################
.-------------------------------------------------------------------.
| Ascii plot of probability distribution for number of chgpts (ncp) |
.-------------------------------------------------------------------.
|Pr(ncp = 0 )=0.009|*                                               |
|Pr(ncp = 1 )=0.991|*********************************************** |
.-------------------------------------------------------------------.
|    Summary for number of Trend ChangePoints (tcp)                 |
.-------------------------------------------------------------------.
|ncp_max    = 1    | MaxTrendKnotNum: A parameter you set           |
|ncp_mode   = 1    | Pr(ncp= 1)=0.99: There is a 99.1% probability  |
|                  | that the trend component has  1 changepoint(s).|
|ncp_mean   = 0.99 | Sum{ncp*Pr(ncp)} for ncp = 0,...,1             |
|ncp_pct10  = 1.00 | 10% percentile for number of changepoints      |
|ncp_median = 1.00 | 50% percentile: Median number of changepoints  |
|ncp_pct90  = 1.00 | 90% percentile for number of changepoints      |
.-------------------------------------------------------------------.
| List of probable trend changepoints ranked by probability of      |
| occurrence: Please combine the ncp reported above to determine    |
| which changepoints below are  practically meaningful              |
'-------------------------------------------------------------------'
|tcp#              |time (cp)                  |prob(cpPr)          |
|------------------|---------------------------|--------------------|
|1                 |3748.000000                |0.84592             |
.-------------------------------------------------------------------.

Of the two possible signal structures--no changepoint(i.e., number of changepoint (ncp)=0 ) vs 1 changepoint (ncp=1), there is strong evidence supporting ncp=1 (the posterior probability is 0.991). You can also plot it out using
 barplot(o$trend$ncpPr, names.arg = c('numer of changepoints=0','numer of changepoints=1'))


