# compute true values monte carlo simulation

#### alex76

##### New Member
Hello,

I'm reading an article in which a monte carlo simulation is performed, but have the following problem:

the data is generated as:
x1=N(2,2)
x2=U(3,8)
P(A=1)=.5
P(Y=1|A,x1,x2)=1/1+exp[-(1.2A-5*x1^2+2*x2)]

the true values are given as : P(Y1=1)=.372 and P(Y0=1)=.352
(where Y1 is the value of Y given that A=1 and Y0 likewise when A=0)

I can't figure out how to compute these true values, so if somebody can give me some help here?

#### Dason

Are you asking how to mathematically figure out what they are or you asking how to do monte carlo simulation to figure out what they are?

#### alex76

##### New Member
I am asking how to figure them out mathematically.
I initially thought: the mean of x1 and x2 are 2 and 5.5 respectively,
so for example P(Y1=1)=1/1+exp[-(1.2-5*2^2+2*5.5) but that's not the idea I guess?

#### Dason

Nope. I'm guessing x1 and x2 and A are all independent?

By the way: I don't think this will be an easy problem? Using monte carlo simulation will get you those answers (or at least something really close to the answers) pretty easily.

#### alex76

##### New Member
Yes, I think they are independent, at least no covariance matrix is given..

I see what you mean by approaching the true values with a monte carlo simulation, but the idea here is that the article generates data according to the laws I just described and discusses several ways to estimate P(Y1=1) and P(Y0=1) or measures following from it (e.g. RR, odds ratio etc)..

So to know how well each method approaches the true values you need to know the actual true values, but how? I need to know since I want to carry out a MC simulation myself with the same goal (comparing different methods)

#### BGM

##### TS Contributor
You are given
$$X_1 \sim N(2, 2)$$
$$X_2 \sim U(3, 8)$$
$$Pr\{Y = 1|A = a, X_1 = x_1, X_2 = x_2\} = \frac {1} {1 + exp\{-(1.2a - 5x_1^2 + 2x_2)\}}$$

If you want, you can convert to a numerical integration problem like this:

$$Pr\{Y = 1|A = a\}$$

$$= \int_{-\infty}^{+\infty} \int_3^8 Pr\{Y = 1|A = a, X_1 = x_1, X_2 = x_2\} f_{X_1}(x_1)f_{X_2}(x_2) dx_2dx_1$$

$$= \frac {1} {5} \int_{-\infty}^{+\infty} \frac {1} {\sqrt{2\pi\sigma^2}} exp\left\{\frac {-(x_1 - 2)^2} {2\sigma^2}\right\} \int_3^8 \frac {1} {1 + exp\{-(1.2a - 5x_1^2 + 2x_2)\}} dx_2dx_1$$

$$= \frac {1} {10} \int_{-\infty}^{+\infty} \left [ln(1 + exp\{-1.2a + 5x_1^2 + 16\}) - ln(1 + exp\{-1.2a + 5x_1^2 + 6\}) \right]$$ $$\frac {1} {\sqrt{2\pi\sigma^2}} exp\left\{\frac {-(x_1 - 2)^2} {2\sigma^2}\right\}dx_1$$

#### Dragan

##### Super Moderator
I am asking how to figure them out mathematically.
I initially thought: the mean of x1 and x2 are 2 and 5.5 respectively,
so for example P(Y1=1)=1/1+exp[-(1.2-5*2^2+2*5.5) but that's not the idea I guess?

No, your idea is not the correct way to do this...Rather, here is how to do this, using Monte Carlo techniques...which is what I think this is really all about.

(1) Generate pseudo-random deviates for X1 and X2 with the specified parameters as given above.

(2) Let Z1 = -(1.2 - 5*X1^2 + 2*X2) for Y = 1.

(3) Let Z2 = -( -5*X1^2 + 2*X2) for Y = 0.

(4) Use the Logistic CDF (Mean 0 and Scale 1) and obtain the the cumulative probabilities for both Z1 and Z2.

(5) Separately compute the averages of the cumulative prob's associated with Z1 and Z2 and these will be .628 and .648, respectively.

(6) Take the complements of the averages (i.e. 1 - .628; and 1 - .648) you compute in Step (5) and these values will yield the answers you are looking for .372 and .352.

And, that should do it.

#### alex76

##### New Member
Thanks again for your advice Dragan, that was indeed what I was looking for. I was wondering, suppose that A is a function of x1 and x2 as well, so that for example now we have:

x1=N(2,2)
x2=U(3,8)

P(A=1|x1,x2)=1/1+exp[-(-(2*x1+2*x2)]

P(Y=1|A,x1,x2)=1/1+exp[-(1.2A-5*x1^2+2*x2)]

Would this change a lot in the procedure?

#### Dragan

##### Super Moderator
Thanks again for your advice Dragan, that was indeed what I was looking for. I was wondering, suppose that A is a function of x1 and x2 as well, so that for example now we have:

x1=N(2,2)
x2=U(3,8)

P(A=1|x1,x2)=1/1+exp[-(-(2*x1+2*x2)]

P(Y=1|A,x1,x2)=1/1+exp[-(1.2A-5*x1^2+2*x2)]

Would this change a lot in the procedure?

I'll think about it, but, before I do, can you make sure to write the correct form of the first expression:

P(A=1|x1,x2)=1/1+exp[-(-(2*x1+2*x2)]

You've got unbalanced parentheses and why do you have two negative signs??

#### Dason

Well if A can only be 0 or 1 why don't you just generate a random observation for A and then for x1 and x2, and then use that to generate an observation to see if Y is 0 or 1. Do that a couple hundred thousand times and compute the average of the Y's. There you have an estimate of the probability Y=1.

#### alex76

##### New Member
Hi Dragan,

yes sorry, the expression should be:

P(A=1|x1,x2)=1/1+exp[-(2*x1+2*x2)]

regards Alex

#### alex76

##### New Member
I think I found the solution. I did it with R, for example the previous example with P(A=1)=.5 also delivered the right values with:

n=10000
n.sim=5000
y1=rep(NA,n.sim)
y0=rep(NA,n.sim)
for (s in 1:n.sim){
x1=rnorm(n,2,2)
x2=runif(n,3,8)
A=rbinom(n,1,.5)
y=plogis(1.2*A-5*x1^2+2*x2)
dat=as.data.frame(cbind(A,y))
y1=mean(dat$y[A==1]) y0=mean(dat$y[A==0])
}
mean(y1)
mean(y0)

so I guess changing A=rbinom(n,1,.5) to A=rbinom(n,1,plogis(2*x1 + 2*x2)) does the trick right?

#### Dragan

##### Super Moderator
I think I found the solution. I did it with R, for example the previous example with P(A=1)=.5 also delivered the right values with:

n=10000
n.sim=5000
y1=rep(NA,n.sim)
y0=rep(NA,n.sim)
for (s in 1:n.sim){
x1=rnorm(n,2,2)
x2=runif(n,3,8)
A=rbinom(n,1,.5)
y=plogis(1.2*A-5*x1^2+2*x2)
dat=as.data.frame(cbind(A,y))
y1=mean(dat$y[A==1]) y0=mean(dat$y[A==0])
}
mean(y1)
mean(y0)

so I guess changing A=rbinom(n,1,.5) to A=rbinom(n,1,plogis(2*x1 + 2*x2)) does the trick right?

Yeah, I think that looks right. What did you get for solutions?...I'll try it on my own using another software package and see how things look.

#### alex76

##### New Member
I got these values:

> mean(y1)
[1] 0.3722326
> mean(y0)
[1] 0.05885827

#### Dragan

##### Super Moderator
I got these values:

> mean(y1)
[1] 0.3722326
> mean(y0)
[1] 0.05885827

I believe you have it correct. I'm still working on it, but, your solution for mean(y1)=.3722326 seems to me to be very reasonable.

#### alex76

##### New Member
Hi there,

In continuation to what I posted before in this thread, I'm still confused about the following:

suppose data is generated as:

W1=rnorm(n,0,1)
W2=runif(n,-2,2)
W3=rnorm(n,1,1)
A <- rbinom(n,1, .5)
Y <- 3.5*A - .75*W1 - 1.5*W2 - .75*W3 + 1.75*(A*W3) + rnorm(n)

where A is a treatment variable. Then the true treatment effect is 5.25, which can be seen straightaway from Y (3.5+1.75) or can be computed as shown in previous posts by MC simulation.

Now, when A is changed to :
A <- rbinom(n,1, 1/(1+exp(-(.4*W1 + .5*W2 + .4*W3))))
it is not as straightforward right?

so I figured the treatment effect can be computed as follows (in R):

n=10000
n.sim=500
y1 =rep(NA,n.sim)
y0 =rep(NA,n.sim)
for (s in 1:n.sim){
W1=rnorm(n,0,1)
W2=runif(n,-2,2)
W3=rnorm(n,1,1)
A <- rbinom(n,1, 1/(1+exp(-(.4*W1 + .5*W2 + .4*W3))))
Y <- 3.5*A - .75*W1 - 1.5*W2 - .75*W3 + 1.75*(A*W3) + rnorm(n)
dat=as.data.frame(cbind(A,Y))
y1=mean(dat$Y[A==1]) y0=mean(dat$Y[A==0])
}
mean(y1)-mean(y0)

which gives a treatment effect of approximately 4.08

however, I'm not sure this is correct??