0
n=50
p=0.32
P=matrix( c(p, 1-p, 0, 0, 0, 0,
p, 0, 1-p, 0, 0, 0,
p, 0, 0, 1-p, 0, 0,
0, p, 0, 0, 1-p, 0,
0, 0, p, 0, 0, 1-p,
0, 0, 0, p, 0, 1-p),
ncol=6, nrow=6, byrow = T)    

X=2  
for(j in 1:n)
{Y=runif(1)      
k=P[X[j],]      
k=cumsum(k)   
if(Y<=k[1])     
{X[j+1]=1}
else if (Y<=k[2])
{X[j+1]=2}
else if (Y<=k[3])
{X[j+1]=3}
else if (Y<=k[4])
{X[j+1]=4}
else if (Y<=k[5])
{X[j+1]=5}
else {X[j+1]=6}}

mean(X)

x=c(1,2,3,4,5,6)
y=c(0.1,0.15,0.22,0.29,0.38,0.45)
approx(x,y,xout=mean(X))

I used the code above to get a mean(y) by linear interpolation with a fixed p. But now, how to change the code in order to plot a graph of mean(y) against p[0:1]??? I kept getting only one mean(y),help me please. p.s.I need only approx$y, that's where I'm stuck :(

Ys Kee
  • 145
  • 1
  • 5

1 Answers1

0

I think it would be better to change your code into function() and use it with sapply(p.vector, ...).

the function
func <- function(p) {
  P = matrix( c(p, 1-p, 0, 0, 0, 0,
                p, 0, 1-p, 0, 0, 0,
                p, 0, 0, 1-p, 0, 0,
                0, p, 0, 0, 1-p, 0,
                0, 0, p, 0, 0, 1-p,
                0, 0, 0, p, 0, 1-p),
              ncol=6, nrow=6, byrow = T)
  X = 2
  for(j in 1:n)
    {Y=runif(1)      
    k=P[X[j],]      
    k=cumsum(k)   
    if(Y<=k[1])     
    {X[j+1]=1}
    else if (Y<=k[2])
    {X[j+1]=2}
    else if (Y<=k[3])
    {X[j+1]=3}
    else if (Y<=k[4])
    {X[j+1]=4}
    else if (Y<=k[5])
    {X[j+1]=5}
    else {X[j+1]=6}}
  return(approx(x, y, xout = mean(X))$y)
}
use the function with p[0:1]
p.vec <- seq(0, 1, 0.01)                     # preparation of p[0:1] as a vector
n = 50                                       # defining other paramaters
x = c(1, 2, 3, 4, 5, 6)
y = c(0.1, 0.15, 0.22, 0.29, 0.38, 0.45)

y.vec <- sapply(p.vec, func)                 # calculation of the y about p[0:1]

plot(p.vec, y.vec, type="o")   # for example

enter image description here

cuttlefish44
  • 6,586
  • 2
  • 17
  • 34