+ Reply to Thread
Page 1 of 2 1 2 LastLast
Results 1 to 15 of 17

Thread: compute true values monte carlo simulation

  1. #1
    Points: 1,647, Level: 23
    Level completed: 47%, Points required for next Level: 53

    Posts
    14
    Thanks
    0
    Thanked 0 Times in 0 Posts

    compute true values monte carlo simulation




    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?

  2. #2
    Devorador de queso
    Points: 95,940, Level: 100
    Level completed: 0%, Points required for next Level: 0
    Awards:
    Posting AwardCommunity AwardDiscussion EnderFrequent Poster
    Dason's Avatar
    Location
    Tampa, FL
    Posts
    12,937
    Thanks
    307
    Thanked 2,630 Times in 2,246 Posts
    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?

  3. #3
    Points: 1,647, Level: 23
    Level completed: 47%, Points required for next Level: 53

    Posts
    14
    Thanks
    0
    Thanked 0 Times in 0 Posts
    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?

  4. #4
    Devorador de queso
    Points: 95,940, Level: 100
    Level completed: 0%, Points required for next Level: 0
    Awards:
    Posting AwardCommunity AwardDiscussion EnderFrequent Poster
    Dason's Avatar
    Location
    Tampa, FL
    Posts
    12,937
    Thanks
    307
    Thanked 2,630 Times in 2,246 Posts
    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.

  5. #5
    Points: 1,647, Level: 23
    Level completed: 47%, Points required for next Level: 53

    Posts
    14
    Thanks
    0
    Thanked 0 Times in 0 Posts
    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)

  6. #6
    TS Contributor
    Points: 22,410, Level: 93
    Level completed: 6%, Points required for next Level: 940

    Posts
    3,020
    Thanks
    12
    Thanked 565 Times in 537 Posts
    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

  7. #7
    Super Moderator
    Points: 13,151, Level: 74
    Level completed: 76%, Points required for next Level: 99
    Dragan's Avatar
    Location
    Illinois, US
    Posts
    2,014
    Thanks
    0
    Thanked 223 Times in 192 Posts
    Quote Originally Posted by alex76 View Post
    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.

  8. #8
    Points: 1,647, Level: 23
    Level completed: 47%, Points required for next Level: 53

    Posts
    14
    Thanks
    0
    Thanked 0 Times in 0 Posts
    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?

  9. #9
    Super Moderator
    Points: 13,151, Level: 74
    Level completed: 76%, Points required for next Level: 99
    Dragan's Avatar
    Location
    Illinois, US
    Posts
    2,014
    Thanks
    0
    Thanked 223 Times in 192 Posts
    Quote Originally Posted by alex76 View Post
    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??

  10. #10
    Devorador de queso
    Points: 95,940, Level: 100
    Level completed: 0%, Points required for next Level: 0
    Awards:
    Posting AwardCommunity AwardDiscussion EnderFrequent Poster
    Dason's Avatar
    Location
    Tampa, FL
    Posts
    12,937
    Thanks
    307
    Thanked 2,630 Times in 2,246 Posts
    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.

  11. #11
    Points: 1,647, Level: 23
    Level completed: 47%, Points required for next Level: 53

    Posts
    14
    Thanks
    0
    Thanked 0 Times in 0 Posts
    Hi Dragan,

    yes sorry, the expression should be:

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

    regards Alex

  12. #12
    Points: 1,647, Level: 23
    Level completed: 47%, Points required for next Level: 53

    Posts
    14
    Thanks
    0
    Thanked 0 Times in 0 Posts
    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[s]=mean(dat$y[A==1])
    y0[s]=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?

  13. #13
    Super Moderator
    Points: 13,151, Level: 74
    Level completed: 76%, Points required for next Level: 99
    Dragan's Avatar
    Location
    Illinois, US
    Posts
    2,014
    Thanks
    0
    Thanked 223 Times in 192 Posts
    Quote Originally Posted by alex76 View Post
    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[s]=mean(dat$y[A==1])
    y0[s]=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.

  14. #14
    Points: 1,647, Level: 23
    Level completed: 47%, Points required for next Level: 53

    Posts
    14
    Thanks
    0
    Thanked 0 Times in 0 Posts
    I got these values:

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

  15. #15
    Super Moderator
    Points: 13,151, Level: 74
    Level completed: 76%, Points required for next Level: 99
    Dragan's Avatar
    Location
    Illinois, US
    Posts
    2,014
    Thanks
    0
    Thanked 223 Times in 192 Posts

    Quote Originally Posted by alex76 View Post
    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.

+ Reply to Thread
Page 1 of 2 1 2 LastLast

           




Similar Threads

  1. monte carlo simulation
    By mathsohard in forum R
    Replies: 1
    Last Post: 03-05-2011, 09:19 PM
  2. Replies: 2
    Last Post: 02-15-2011, 06:37 AM
  3. Replies: 2
    Last Post: 03-21-2010, 01:25 AM
  4. What is a monte carlo simulation?
    By brad henrie in forum Statistical Research
    Replies: 1
    Last Post: 10-22-2009, 08:53 PM
  5. Replies: 4
    Last Post: 07-23-2009, 02:03 AM

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts






Advertise on Talk Stats