# BS1033 Lecture 4 Analysis
# Author: Chris Hansman
# Email: [email protected]
# Date : 28/01/19
# Installing Packages
#install.packages(tidyverse)
# Loading Libraries
library(tidyverse)
##
# Functions in r
##
# A function that adds two numbers
addition <- function(x,y){x+y}#Evaluating Our Addition Functionaddition(2,2)# A function that squares a number minus 1xsquared <- function(x){(x-1)^2}#Evaluating our Square Functionxsquared(1)# Specify range of x-axisggplot(data.frame(x = c(0,3) ), aes(x=x)) + stat_function(fun = xsquared) #Minimizing our square functionz<-optim(par=1, fn=xsquared)z$par#————————————————–## Simulating Probit Data#————————————————–##Setting Seedset.seed(050187)# Choosing parametersbeta_0 <- 0.2beta_1 <- 0.5n <- 10000# Simulating Datadata <- tibble(#Simulating x and vx_i = rnorm(n),v_i = rnorm(n),#Generating latent indexy_i_star=beta_0+beta_1*x_i+v_i,#Generating observed binary variabley_i = y_i_star>0
)
#Converting Data to Vector Form
x_i <- data.matrix(data$x_i)y_i <- data.matrix(data$y_i)#————————————————–## Writing a Likelihood Function#————————————————–#loglik <- function(beta){#Log(Phi(Xbeta))l_1 <- log(pnorm(beta[1]+beta[2]*x_i))#Log(1-Phi(Xbeta))l_2 <- log(1-pnorm(beta[1]+beta[2]*x_i))#Summing Everything-sum(y_i*l_1+(1-y_i)*l_2)}loglik(c(1,1))#Minimizing our likelihood functionprobit_ll<-optim(par=c(1,1), fn=loglik)#Value of Likelihood at minimumprobit_ll$val#Value of Parameters at Minimumprobit_ll$par# The Easy Wayprobit_glm <- glm(y_i ~ x_i, family = binomial(link = “probit”), data = data)summary(probit_glm)logLik(probit_glm)#Plotting Likelihood#Creating Univariate Functionloglik_b1 <- function(x){(loglik(c(0.2,x)))} #Vectorizing Univariate Functionloglik_b1<-Vectorize(loglik_b1, vectorize.args=”x”)#Vectorizing Univariate Function#Plotting Likelihood#Creating Univariate Functionloglik_b2 <- function(x){(-loglik(c(0.2,x)))} #Vectorizing Univariate Functionloglik_b2<-Vectorize(loglik_b2, vectorize.args=”x”)#Plotting Likelihoodggplot(data.frame(x = c(0.25,0.75) ), aes(x=x)) + stat_function(fun = loglik_b1) +ylab(“Log-Likelihood”) +xlab(expression(beta[1])) +theme_classic()ggsave(“Log_Likelihood.pdf”)#Plotting Likelihoodggplot(data.frame(x = c(0.25,0.75) ), aes(x=x)) + stat_function(fun = loglik_b2) +ylab(“Log-Likelihood”) +xlab(expression(beta[1])) +theme_classic()ggsave(“Log_Likelihood2.pdf”)#Plotting Likelihoodggplot(data = data.frame(x = c(-5, 5)), aes(x)) +stat_function(fun = pnorm, n = 101, args = list(mean = 0, sd = 1)) + ylab(“Probability”) +theme_classic()ggsave(“cdf.pdf”)#Plotting Likelihoodggplot(data = data.frame(x = c(-5, 5)), aes(x)) +stat_function(fun = plogis, n = 101) + ylab(“Probability”) +theme_classic()ggsave(“cdf2.pdf”)
Reviews
There are no reviews yet.