## Loading required package: carData
## Loading required package: sp
## Loading required package: raster
##
## Attaching package: 'raster'
## The following objects are masked from 'package:MASS':
##
## area, select
## Checking rgeos availability: TRUE
## rgdal: version: 1.5-16, (SVN revision 1050)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: C:/Users/lfult/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: C:/Users/lfult/Documents/R/win-library/4.0/rgdal/proj
## Linking to sp version:1.4-2
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
## rgeos version: 0.5-5, (SVN revision 640)
## GEOS runtime version: 3.8.0-CAPI-1.13.1
## Linking to sp version: 1.4-2
## Polygon checking: TRUE
##
## Attaching package: 'shiny'
## The following object is masked from 'package:ggExtra':
##
## runExample
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
## Loading required package: Matrix
## Registered S3 methods overwritten by 'spatialreg':
## method from
## residuals.stsls spdep
## deviance.stsls spdep
## coef.stsls spdep
## print.stsls spdep
## summary.stsls spdep
## print.summary.stsls spdep
## residuals.gmsar spdep
## deviance.gmsar spdep
## coef.gmsar spdep
## fitted.gmsar spdep
## print.gmsar spdep
## summary.gmsar spdep
## print.summary.gmsar spdep
## print.lagmess spdep
## summary.lagmess spdep
## print.summary.lagmess spdep
## residuals.lagmess spdep
## deviance.lagmess spdep
## coef.lagmess spdep
## fitted.lagmess spdep
## logLik.lagmess spdep
## fitted.SFResult spdep
## print.SFResult spdep
## fitted.ME_res spdep
## print.ME_res spdep
## print.lagImpact spdep
## plot.lagImpact spdep
## summary.lagImpact spdep
## HPDinterval.lagImpact spdep
## print.summary.lagImpact spdep
## print.sarlm spdep
## summary.sarlm spdep
## residuals.sarlm spdep
## deviance.sarlm spdep
## coef.sarlm spdep
## vcov.sarlm spdep
## fitted.sarlm spdep
## logLik.sarlm spdep
## anova.sarlm spdep
## predict.sarlm spdep
## print.summary.sarlm spdep
## print.sarlm.pred spdep
## as.data.frame.sarlm.pred spdep
## residuals.spautolm spdep
## deviance.spautolm spdep
## coef.spautolm spdep
## fitted.spautolm spdep
## print.spautolm spdep
## summary.spautolm spdep
## logLik.spautolm spdep
## print.summary.spautolm spdep
## print.WXImpact spdep
## summary.WXImpact spdep
## print.summary.WXImpact spdep
## predict.SLX spdep
##
## Attaching package: 'spdep'
## The following objects are masked from 'package:spatialreg':
##
## anova.sarlm, as.spam.listw, as_dgRMatrix_listw, as_dsCMatrix_I,
## as_dsCMatrix_IrW, as_dsTMatrix_listw, bptest.sarlm, can.be.simmed,
## cheb_setup, coef.gmsar, coef.sarlm, coef.spautolm, coef.stsls,
## create_WX, deviance.gmsar, deviance.sarlm, deviance.spautolm,
## deviance.stsls, do_ldet, eigen_pre_setup, eigen_setup, eigenw,
## errorsarlm, fitted.gmsar, fitted.ME_res, fitted.sarlm,
## fitted.SFResult, fitted.spautolm, get.ClusterOption,
## get.coresOption, get.mcOption, get.VerboseOption,
## get.ZeroPolicyOption, GMargminImage, GMerrorsar, griffith_sone,
## gstsls, Hausman.test, HPDinterval.lagImpact, impacts, intImpacts,
## Jacobian_W, jacobianSetup, l_max, lagmess, lagsarlm, lextrB,
## lextrS, lextrW, lmSLX, logLik.sarlm, logLik.spautolm, LR.sarlm,
## LR1.sarlm, LR1.spautolm, LU_prepermutate_setup, LU_setup,
## Matrix_J_setup, Matrix_setup, mcdet_setup, MCMCsamp, ME, mom_calc,
## mom_calc_int2, moments_setup, powerWeights, predict.sarlm,
## predict.SLX, print.gmsar, print.ME_res, print.sarlm,
## print.sarlm.pred, print.SFResult, print.spautolm, print.stsls,
## print.summary.gmsar, print.summary.sarlm, print.summary.spautolm,
## print.summary.stsls, residuals.gmsar, residuals.sarlm,
## residuals.spautolm, residuals.stsls, sacsarlm, SE_classic_setup,
## SE_interp_setup, SE_whichMin_setup, set.ClusterOption,
## set.coresOption, set.mcOption, set.VerboseOption,
## set.ZeroPolicyOption, similar.listw, spam_setup, spam_update_setup,
## SpatialFiltering, spautolm, spBreg_err, spBreg_lag, spBreg_sac,
## stsls, subgraph_eigenw, summary.gmsar, summary.sarlm,
## summary.spautolm, summary.stsls, trW, vcov.sarlm, Wald1.sarlm
## The following objects are masked from 'package:elsa':
##
## geary, moran
## -- Attaching packages --------------------------------------------- tidyverse 1.3.0 --
## v tibble 3.0.3 v dplyr 1.0.2
## v tidyr 1.1.2 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.5.0
## v purrr 0.3.4
## -- Conflicts ------------------------------------------------ tidyverse_conflicts() --
## x dplyr::combine() masks gridExtra::combine()
## x tidyr::expand() masks Matrix::expand()
## x tidyr::extract() masks raster::extract()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x tidyr::pack() masks Matrix::pack()
## x dplyr::recode() masks car::recode()
## x dplyr::select() masks raster::select(), MASS::select()
## x purrr::some() masks car::some()
## x tidyr::unpack() masks Matrix::unpack()
## Loading required package: lattice
##
## Attaching package: 'caret'
## The following object is masked from 'package:purrr':
##
## lift
## Loaded glmnet 4.0-2
##
## Attaching package: 'spatialEco'
## The following object is masked from 'package:dplyr':
##
## combine
## The following object is masked from 'package:spData':
##
## elev
## The following object is masked from 'package:gridExtra':
##
## combine
## The following object is masked from 'package:grid':
##
## explode
## The following object is masked from 'package:elsa':
##
## correlogram
## The following object is masked from 'package:raster':
##
## shift
## Registered S3 method overwritten by 'gplots':
## method from
## reorder.factor gdata
## Registered S3 methods overwritten by 'lme4':
## method from
## cooks.distance.influence.merMod car
## influence.merMod car
## dfbeta.influence.merMod car
## dfbetas.influence.merMod car
## Loaded lars 1.2
## Loading required package: survival
##
## Attaching package: 'survival'
## The following object is masked from 'package:spatialEco':
##
## concordance
## The following object is masked from 'package:caret':
##
## cluster
#Load shape and merge with data
myshape=shapefile("D:/Covid19/cb_2018_us_county_500k")
countydata=read.csv("D:/Covid19/reducedcounty4.csv",fileEncoding="UTF-8-BOM")
countydata$M=as.numeric(countydata$countyFIPS)
myshape$M=as.numeric(myshape$GEOID)
counties=merge(myshape, countydata, by="M",type="left")
#deaths=read.csv("https://usafactsstatic.blob.core.windows.net/public/data/covid-19/covid_deaths_usafacts.csv", fileEncoding="UTF-8-BOM")
#deaths$Mean=rowMeans(deaths[, 5:ncol(deaths)])
#deaths$Var=apply(deaths[, 5:ncol(deaths)],1, var)
#mydata=deaths[deaths$Var<deaths$Mean,]
#mydata$M=as.numeric(mydata$countyFIPS)
#write.csv(mydata, "c:/users/lfult/desktop/ram/deaths.csv", row.names = FALSE)
mydata=read.csv("c:/users/lfult/desktop/ram/deaths.csv")
mydata$theta1=round(mydata$Mean^(3/2)/mydata$Var,3)
mydata$gamma1=round(abs((1-sqrt(mydata$Mean/mydata$Var))/mydata$theta1),3)
#Moment Estimator Plots of Theta and 1 / Gamma
plot(density(mydata$theta1))
plot(density(1/mydata$gamma1), xlim=c(0,100))
#p-value
mydata$value=rep(0, nrow(mydata))
for (i in 1:nrow(mydata)){
for(j in 9:190){p1=
sum((mydata[i,j]-1)*log(1+mydata$gamma1[i]*mydata[i,j]))
}
mydata$value[i]=
-2*mydata$theta1[i]*mydata$gamma1[i]*mydata$Mean[i]+
182*(mydata$theta1[i]-mydata$Mean[i])-
182*mydata$Mean[i]*log(mydata$theta1[i]/mydata$gamma1[i])-
p1
}
mydata$p=round(1-pchisq(mydata$value,1),3)
notsig=mydata$County.Name[mydata$p>.05]
notsig2=mydata$State[mydata$p>.05]
paste(notsig,notsig2)
## [1] "Kenai Peninsula Borough AK" "Matanuska-Susitna Borough AK"
## [3] "Conway County AR" "Lafayette County AR"
## [5] "Logan County AR" "Van Buren County AR"
## [7] "Mariposa County CA" "Mono County CA"
## [9] "Delta County CO" "Elbert County CO"
## [11] "Huerfano County CO" "Kit Carson County CO"
## [13] "La Plata County CO" "Pitkin County CO"
## [15] "Saguache County CO" "Teller County CO"
## [17] "Clay County GA" "Schley County GA"
## [19] "Webster County GA" "Gooding County ID"
## [21] "Effingham County IL" "Blackford County IN"
## [23] "Brown County IN" "Fountain County IN"
## [25] "Fulton County IN" "Huntington County IN"
## [27] "Jasper County IN" "Owen County IN"
## [29] "Washington County IN" "Appanoose County IA"
## [31] "Buchanan County IA" "Butler County IA"
## [33] "Clayton County IA" "Madison County IA"
## [35] "Taylor County IA" "Van Buren County IA"
## [37] "Bourbon County KS" "Clay County KS"
## [39] "Dickinson County KS" "Jackson County KS"
## [41] "Kearny County KS" "Meade County KS"
## [43] "Breckinridge County KY" "Henry County KY"
## [45] "Knott County KY" "McLean County KY"
## [47] "Meade County KY" "Metcalfe County KY"
## [49] "Pike County KY" "Nantucket County MA"
## [51] "Alcona County MI" "Cheboygan County MI"
## [53] "Dickinson County MI" "Emmet County MI"
## [55] "Gladwin County MI" "Gogebic County MI"
## [57] "Mecosta County MI" "Brown County MN"
## [59] "Chippewa County MN" "Chisago County MN"
## [61] "Kandiyohi County MN" "Todd County MN"
## [63] "Wilkin County MN" "Franklin County MS"
## [65] "Bates County MO" "Callaway County MO"
## [67] "Carter County MO" "Lewis County MO"
## [69] "Lincoln County MO" "Shannon County MO"
## [71] "Statewide Unallocated MT" "Fergus County MT"
## [73] "Glacier County MT" "Lincoln County MT"
## [75] "Sweet Grass County MT" "Dixon County NE"
## [77] "Furnas County NE" "Lincoln County NE"
## [79] "Polk County NE" "Richardson County NE"
## [81] "Saline County NE" "Lander County NV"
## [83] "Cheshire County NH" "Quay County NM"
## [85] "Allegany County NY" "Cayuga County NY"
## [87] "Oswego County NY" "Schoharie County NY"
## [89] "Wayne County NY" "Dare County NC"
## [91] "Perquimans County NC" "Ramsey County ND"
## [93] "Walsh County ND" "Athens County OH"
## [95] "Brown County OH" "Champaign County OH"
## [97] "Morrow County OH" "Union County OH"
## [99] "Choctaw County OK" "Cotton County OK"
## [101] "Craig County OK" "Kiowa County OK"
## [103] "Latimer County OK" "Logan County OK"
## [105] "Major County OK" "Pawnee County OK"
## [107] "Pontotoc County OK" "Roger Mills County OK"
## [109] "Tillman County OK" "Josephine County OR"
## [111] "Clarion County PA" "Fulton County PA"
## [113] "McKean County PA" "Mifflin County PA"
## [115] "Tioga County PA" "Lewis County TN"
## [117] "Lincoln County TN" "Moore County TN"
## [119] "Cottle County TX" "Crane County TX"
## [121] "Mitchell County TX" "Iron County UT"
## [123] "Addison County VT" "Washington County VT"
## [125] "Windham County VT" "Windsor County VT"
## [127] "Gloucester County VA" "Greene County VA"
## [129] "New Kent County VA" "Klickitat County WA"
## [131] "Lincoln County WA" "Stevens County WA"
## [133] "Clay County WV" "Lewis County WV"
## [135] "Marion County WV" "Roane County WV"
## [137] "Buffalo County WI" "Columbia County WI"
## [139] "Juneau County WI" "Kewaunee County WI"
## [141] "Marquette County WI" "Monroe County WI"
## [143] "Polk County WI" "Sauk County WI"
## [145] "Campbell County WY" "Carbon County WY"
table(notsig2)
## notsig2
## AK AR CA CO GA IA ID IL IN KS KY MA MI MN MO MS MT NC ND NE NH NM NV NY OH OK
## 2 4 2 8 3 7 1 1 8 6 7 1 7 6 6 1 5 2 2 6 1 1 1 5 5 11
## OR PA TN TX UT VA VT WA WI WV WY
## 1 5 3 3 1 3 4 3 8 4 2
counties1=merge(counties,mydata, by="M",type="left")
counties2=sp.na.omit((counties1))
## Deleting rows: 2728293032232332432532632732848348450250350450550650750850951051175284688288395710681069107612001201120212031210139614301431144214431475147614771487151515261527152815841592160316041615161616171618161916251649172617331777186018641882188519091910191319321933195420222041205420622063206820822112211721692182218521862258242725392541257725782579258825922623263526592660266727092716272527432748276327862792280828442890289128972899290129343054345689101112131415161718192021222324252631323334353637383940414243444546474849505154555758596061626364656667686970717273747576777879808182838485868788899091939495969899100101102104105106107108109110111112114115116117118120121122123124129130131132133134135136137138139140141142145146147148150151152153154155156157158160162163164166167168169170171172173174175176178179180181182183184185186187188189190191192194195196198200201202203204205206207208209210211212213214215216217218219220222223224225226227228229230231233234235236237238239240241243244245246247248249250251252253254255257258259260261262263264265266267268269270271272273274276277278279280281282283284285286287288289290291292293295296297299300302303304305306307308309310311312313314316317318319321329330331332333335336337338340341342343344345346347349350352354355356357358359360361362363364365366367368369371372373375376377378381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413415418419420421422424426427428429430431433434435437438440441442443444445446447448449450453454455456457458459460461462463464465466469470472473474475476477478479480481482485486487489490491492493494495497498499501512513514515516517518519520521522524525526527528529530531532533535536537538539540541542543544545546547548549550551552553554555556558559560561562563564565566567568569570572573574575576577578579580581582583584585586587589590592593594595596597598599600601602603604605607608609610611612613614615616617618619620622623625626627628630631632633634635636637638639640641642643644646647648649650651652653655658659660661662663664665666667668669671672673675676677678680681682683684685686687689691692694695696698699700702703704706707709710711712713714715716717718719720721722723724725726727728729730731732733734735736738739740741742743745746747748750751753755757758759760761762763764765766767768769770771772773774775778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817819820821822823824825826827828830831832833834835836837838839840841842845847848850851852854855856857858859860862863864865866867868869870873874875876877878879880881884885886887888889891893894895896897900901902903904905906907908909912914916917919920922923924925926927928929930931932933935936937938939940942943944945946947948950952953954955956958959960961962963964965966967968969970971973974975976978979981982984985987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101310151016101710181019102110221024102510261027102810291030103210331034103510361037103810391040104210431044104510461047104810491050105210531054105510561058105910601061106210631064106610701071107210731074107510771078107910801081108210831085108710891090109110921093109410951096109710981099110011011103110411051106110711081109111011111112111311141115111711181121112211231124112511271128112911301131113311341135113611391140114111421143114411451147114811491152115311541155115811591161116211631164116511661167116811691170117111721174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912041205120612071208120912111212121312141215121612171218121912201221122212231224122512261228122912301231123212331234123512361237123812391240124112421243124412451246124712491250125112521253125512561258125912601262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161318131913201321132213231324132513261327132913301331133313341335133613371338133913401341134213431344134613471348134913501351135213531354135613571358136013611362136313641365136713681369137113721374137513781379138013811382138313841385138713881389139013911392139313941395139713991400140314041405140614071409141014111412141314141415141614171418141914201421142214231424142514271428142914321433143414351436143714391440144414461448144914501451145214541455145614571460146114621463146414651466146714681469147114721473147414781479148114821483148414851488148914901491149214931494149514961497149814991500150215041507150815091510151115121513151615171518151915201522152315241525153015311532153315341536153715381540154115421543154415451546154715481549155115521553155415551556155715581561156215631564156615671568156915711572157315741575157615771581158215831585158715901591159315941595159615971598160016011602160516061607160816091610161116121613161416201621162216231624162716291630163116321633163416351636163716391640164116421643164416451646164716481650165116521653165416551656165716581659166016611662166316651666166716681669167016711672167416751676167716781679168016811682168316841685168616871688168916901692169416961698169917001702170317041705170617071708170917101712171317141715171617171718171917201721172217231724172517271728172917301731173417351736173717381739174017411742174317441745174617471748175017531754175517561758175917601761176217631765176717681769177017711772177317741775177617781779178017811782178317841785178617871788178917901791179217931795179617971798179918011802180318051807180818091810181118121813181418161817181918201821182218231824182518261827182918301831183318341835183618371838183918401841184218431844184518461847184818491850185118521853185418551856185918611862186318651866186718691870187218741875187618771878187918801881188318841887188818891890189318941895189618971898189919011903190419051906190719081912191419151916191719181920192119221923192519261927192819291930193119341935193819391940194119421943194419451946194719481950195119521953195519571958195919611964196519661967196819691970197119721973197519781979198119821983198419851986198719881989199119921993199419951996199719981999200020012002200320042006200720082009201020112012201320142015201620172018202120232024202520262027202820292030203120322033203520362037203820392040204220432044204520462047204820492050205120522056205720592060206120642065206620692070207120722073207420752076207720782079208020812083208420852087208820892090209120922093209420952096209720992100210121022103210421062107210921102111211321142115211821192120212121222123212421252126212721282129213021312134213521362137213821392140214121422143214421452146214721482149215021512152215321542155215721582159216021612162216321642165216621672168217021712172217321742175217621782179218021812183218421872188218921902191219221952196219721982199220122022203220522072208220922102211221222142215221622172218221922212222222322242225222622272228222922312233223422352238223922412242224422452246224722482249225022512253225422552257225922602262226322662267226822702272227322742275227622772278228022812282228322852286228822892290229122922293229422962297229923002301230223032304230523062307230823092310231123122313231423152317231823202321232223242325232623272328232923302332233323352336233723382340234223432345234623472348234923502351235523562357235823592360236123622363236423652366236723682369237023712372237323742375237623772378237923802382238323842385238623872388238923902391239223932394239523962397239823992402240324042406240724082409241124122413241424152416241724192420242124222423242424262428243024312433243424382439244024412442244324442445244624472448245024512452245324542455245624572458246024612463246424662467246824692470247124732474247524762477247824792480248124822483248524862487248824922493249524962497249925002503250425052506250725082509251025112512251325142515251625182519252025212522252325242526252725282529253025322533253425362537253825402544254525472548254925502551255325542555255625572558255925602561256225642565256625672568256925702571257225732575257625802581258225832589259025912593259425952596259725982599260026012602260526062607260826092610261226132614261526162617261826202622262426252626262726282629263026312632263326342636263726382639264026422643264426452646264726482649265126522653265426552656265726582661266226632664266526682670267126722673267426752676267726782679268126822683268426852686268726882689269026922693269426952696269726982699270027012702270327052706270827102711271327142715271827192720272127232724272627272728272927302731273227332734273527362737273827392740274127442745274727512752275327542755275627592760276127622764276527662767276927702771277227732774277527762778277927802781278227832785278827892790279127932794279527972798279928002801280228032804280528062807280928102812281428152816281728182819282028222823282428252826282728282829283028312833283528362837283928402841284228452846284728482849285228532854285628572858285928602861286228632864286528662867286828692870287128722873287428772878287928802881288228832884288528862887288828892892289328942895289628982900290229032904290529072908291029112912291329142915291629172918291929202921292229232924292529262927292829292930293229332935293629372939294029412942294529462948294929502952295429552956295829592960296129622963296429652966296729682969297029722973297429752976297729782980298129822983298429852986298729882989299129922993299429952996299829993000300130023003300430053006300730083009301030113012301330143015301730183019302030223023302430253026302730283029303030313032303330353036303730383039304130423043304430453046304730483049305030513052305330553057305830593061306230643065306630673068306930703071307230733074307530763077308030813082308330843085308630873088308930903091309230933095309630973098309931003101310331043105310631073108310931103111311331143115311631173118311931203121312231233124312531263127312831293130313131333135313731393140314131423143314431463147314831493150315131523153315531563157315831593160316131623163316431653166316831693172317331743175317631773178317931803181318231833184318531863189319031913192319331943195319631973198319932003201320232033204320532063207320832103211321232133214321532163217321832193220322132223224322532263227322832293230323132323233
#############################################################################
library(totalcensus)
counties2@data$ABB=convert_fips_to_names(counties2@data$STATEFP)
qpal<-colorBin("Reds",c(0,1,2,3,4),bins=5,pretty=TRUE, alpha=.1)
leaflet() %>%
#Overlays
addTiles(group = "OSM (default)") %>%
#Panes
addMapPane("borders", zIndex = 410) %>%
#Base Diagrams
addPolylines(data = counties2,
color = "black",
opacity = 1,
weight = 1, group="Borders", options = pathOptions(pane = "borders"))%>%
fitBounds(-124.8, -66.9, 24.4,49.4) %>%
setView(-98.6, 39.83, zoom = 4)%>%
addPolygons(data=counties2,stroke = FALSE,
fillOpacity = 1, smoothFactor = 0.2,
color=~qpal(counties2@data$Mean),
popup = paste(
"County: ", counties@data$NAME, "<br>",
"Theta: ", counties@data$theta1, "<br>",
"Gamma: ", counties2@data$gamma1, "<br>",
"p-value: ", counties2@data$p, "<br>",
"Mean: ", counties2@data$Mean, "<br>",
"Var: ", counties2@data$Var))%>%
addLegend(data=counties2,
"bottomleft", opacity=1, pal = qpal,
values = ~Mean,
title = "Mean Estimates, IRL Poisson")%>%
addLayersControl(
overlayGroups = c("Borders"),
options = layersControlOptions(collapsed = FALSE)
)
#Maximum Likelihood Estimator
#z=rep(0, nrow(mydata))
#subdata=mydata[,8:190]
#n=ncol(subdata)
#m=nrow(subdata)
#library(pracma)
#for (j in 1:m){ #each row has an estimator
# temp=function(theta2) #build a function for Newton Rhapson
# {
# y=as.numeric(subdata[j,]) #select the first row variable
#y_1=y-1 #y-1
# ybar=mean(y) #ybar
# y_m=y-ybar #y-ybar
#yxm=y*ybar #y x ybar
#n=length(y) #n
#sum[(y*y-1)/(abs(theta*(y-ybar))+y*ybar)] - n
#theta2-(1/(n*ybar))*
# sum(abs(y*(y_1))/
# (1+(1/theta2 -1/ybar)*ybar))
#}
#We solve using Simpson's method (Newton Raphson)
#z[j]=suppressWarnings(newton(temp,0 , maxiter = 500, tol = 1e-04)$root)
# }
#plot(density(z))
#From these results, it is clear there is an error in either 1) the programming, 2) the formula, 3) both.