Ralph Probabilidade 1.0

augustofilho — Jun 4, 2014, 6:43 PM

#Construcao do Modulo de Simulacao

banco2=function(){
  #Semente fixa igual a 2014
  set.seed(2014)

  #Carregando pacotes
  require(stats)

  #Tamanho da amostra
  cat("Digite o tamanho da amostra:\n")
  n=scan( ,nlines=1)

  #Entrando com os valores de z - vetor de 0 ate C
  cat("Informe o vetor de z - De 0 ate C")
  z=scan()

  #Simulando o alpha - parametro de discriminação de atitude          
  alpha=seq(0.2,4,length=n) 

  #Simulando o delta - parametro de localização de atitude   
  delta=seq(-2.4,-1.6,length=n) 

  #Simulando theta - parametro da localização do individuo em um continuum latente
  theta=rnorm(n,0,1)               

  #Restricoes para os tau's:
  #artigo do J.S. Roberts, p. 178
  #1) tau0=0  (por definicao)
  #2) tau(z) = tau(-z)   (simetrico em torno do zero)
  #3) soma(tau_{k=0 ate M)=0

  tau=rep(0,3) #cria um vetor com zeros de tamanho 3 para
  #receber os valores tau a esquerda do zero

  tau[1]=runif(1,min=-3,max=-2)
  tau[2]=runif(1,min=-2,max=-1.33)  
  tau[3]=runif(1,min=-1,max=-0.67) 

  #Vetor final de tau's
  tau=c(tau,0,-tau)

  #Agrupando a saida de dados

  cat("Parametros dos itens utilizados na primeira simulacao\n") 
  cat("\n")
  saida=cbind(theta,alpha,delta,tau[1],tau[2],tau[3],tau[4],tau[5],tau[6],tau[7])
  colnames(saida)=c("theta","alpha","delta","tau1","tau2","tau3","tau4","tau5","tau6","tau7")


  cat("Probabilidade *****P(Z=z|theta)****\n") 
  cat("\n")

  #C is the number of observable response categories minus 1

  i=1
  j=1
  for(i in 1:length(theta)){
    for(j in 1:length(theta)){

      C=length(tau)-1
      M=2*C+1.
      NUM_1 = exp(alpha[j]*(z*(theta[i] - delta[j]) - sum(tau[1:(z+1)]))) 
      NUM_2 = exp(alpha[j]*((M - z)*(theta[i] - delta[j]) - sum(tau[1:(z+1)])))
      DEN = 0

      for (w in 0:C){
        DEN=DEN + exp(alpha[j]*(w*(theta[i] - delta[j]) - sum(tau[1:(w+1)]))) + 
          exp(alpha[j]*((M - w)*(theta[i] - delta[j]) - sum(tau[1:(w+1)])))
      }
      resultado=(NUM_1 + NUM_2)/DEN
      cat("\n O resultado é: \n",resultado)
      #return((NUM_1 + NUM_2)/DEN)
    }
  }

}