Scribd
Upload
9 views4 pages

R Functions for Statistical Analysis

This document contains R code for several statistical analysis questions. Question 2 defines a function MAMD that calculates the mean absolute deviation from the median for a vector. Question 3 defines two functions, MNTL2 and MNTL_weak, that simulate time until a company loses power multiple times and return the mean and confidence interval. Question 4 defines a function estimate that takes a sample size, samples from a uniform distribution, and returns the variance of transforming the samples using exponential and power functions to compare their variances as the sample size increases.

Uploaded by

Doritosxu
Copyright
© All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as DOCX, PDF, TXT or read online on Scribd
9 views4 pages

R Functions for Statistical Analysis

This document contains R code for several statistical analysis questions. Question 2 defines a function MAMD that calculates the mean absolute deviation from the median for a vector. Question 3 defines two functions, MNTL2 and MNTL_weak, that simulate time until a company loses power multiple times and return the mean and confidence interval. Question 4 defines a function estimate that takes a sample size, samples from a uniform distribution, and returns the variance of transforming the samples using exponential and power functions to compare their variances as the sample size increases.

Uploaded by

Doritosxu
Copyright
© All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as DOCX, PDF, TXT or read online on Scribd

STAT 340 Code:

Assignment 1 Q2

MAMD=function(v)
{s=NULL #Denote s and output as NULL vectors
output=NULL
q=median(v) #Denote q as the median (2nd quartile) of the
vector
n=length(v) #Denote n as the length (# values) of the vector
for (i in 1:n) # loop n times, find abs(each value - median)
and store in vector s
{s=c(s,abs(v[i]-q))}
output=c(output,sum(s)/n) #divide sum of s vector by n
return (output)}






























Assignment 1 Q3 PA

MNTL2 = function(N) { #N is the number of times we make MNTL lose
power!
output = NULL; #vector which is empty.
interval <- c(0, 0, 0)
# A variable that tells me how many yrs until MNTL lost power.
for (i in 1:N) {#start of for loop
T = rexp(9,1/20)
#9 data values from an exponential distribution with a rate of 1/20
#the r_distribution name is used to create random data, eg.
rnorm, rbinom, runif....
P1 = min(T[1],T[2],T[3],T[4],T[8])
P2 = min(T[1],T[6],T[5],T[4],T[3],T[7])
P3 = min(T[1],T[2],T[7])
P4 = min(T[1],T[6],T[5],T[8])
L = max(P1,P2,P3,P4)
output = c(output,L)
}#End for loop
error=(qnorm(0.975)*sd(output))/sqrt(N)
# denote the standard error as 1.96 * stddev of the N
realizations/sqrt of N
interval[1]<-c(mean(output)-error)
interval[2]<-c(mean(output))
interval[3]<-c(mean(output)+error)
# denote interval as the vector containing the lower limit, mean,
and upper limit of the N realizations with 95% confidence.
print(interval)}
#end function


outputs:

>source("/Users/DavidXu/Documents/Lectures/STAT340/Assignments/A1/A1Q3.
R")
> MNTL2(100)
[1] 7.169840 8.326651 9.483461
>






Assignment 1 Q3 PB

MNTL_weak = function(N) { #N is the number of times we make MNTL
lose power!
output = NULL; #vector which is empty.
interval <- c(0, 0, 0)
# A variable that tells me how many yrs until MNTL lost power.
for (i in 1:N) {#start of for loop
T = rexp(9,1/20)
#9 data values from an exponential distribution with a rate of
1/20
#the r_distribution name is used to create random data, eg.
rnorm, rbinom, runif....
Z=rnorm(6,13,4)
# 6 cities (excluding MNTL itself) that may lose ability to
produce power
P1 = min(T[1],T[2],T[3],T[4],T[8], Z[1], Z[2], Z[4], Z[5], Z[6])
P2 = min(T[1],T[6],T[5],T[4],T[3],T[7], Z[1], Z[2], Z[3], Z[4],
Z[5], Z[6])
P3 = min(T[1],T[2],T[7], Z[1], Z[2], Z[4])
P4 = min(T[1],T[6],T[5],T[8], Z[1], Z[2], Z[3], Z[6])
L = max(P1,P2,P3,P4)
output = c(output,L)
}#End for loop
error=(qnorm(0.975)*sd(output))/sqrt(N)
interval[1]<-c(mean(output)-error)
interval[2]<-c(mean(output))
interval[3]<-c(mean(output)+error)
print(interval)}
#end function

outputs:

>source("/Users/DavidXu/Documents/Lectures/STAT340/Assignments/A1/A1Q/
PB.R")
> MNTL_weak(100)
[1] 5.581896 6.227770 6.873645







Assignment 1 Q4 PA

# By the properties of variance, and given X(i) i.i.d

#V(sum of i, i=1 to n, 2^X(i))= sumof i, i=1 to n V(2^X(i))
#V(sum of i, i=1 to n, exp(X(i)))= sumof i, i=1 to n V(exp(X(i)))
# hence, the variance for both increases as n increases, estimate
V(2^X(i)) and V(exp(X(i))) to compare the variances

estimate=function(n)
{variance=NULL
variance=c(variance, var(2^(runif(n))))
variance=c(variance, var(exp(runif(n))))
print(variance)}

# This function will produce the variance of exp(X(i)) and 2^X(i)
after providing a sample size for X(i) (the larger n is, the more
accurate)


outputs:

>source("/Users/DavidXu/Documents/Lectures/STAT340/Assignments/A1/A1Q4
PA.R")
> estimate(10000000)
[1] 0.08264784 0.24209412

You might also like