Introduction

In the following report, I attempt to answer the following questions regarding crime in Chicago. Given that I work closely with CPD on a daily basis, these data are quite familiar to me.

Question 1)

How has crime changed over time?

Explore the data. Provide some methodology on how you would attack this question. And then provide some Python code which analyzes that data to test your questions or hypotheses.

Question 2)

How has crime changed by type of crime? E.g. how has homicide changed over time? How have other crime types changed over time? What are most remarkable?

Provide some methodology on how you would attack this question. And then provide some Python code which analyzes that data to test your questions or hypotheses.

Question 3)

Can you forecast how crime is likely to change in the future?

Provide some methodology on how you would attack this question. And then provide some Python code which analyzes that data to test your questions or hypotheses.

How has crime changed over time?

When examining how crime has changed over time, there are two main types of analysis that are worth considering.

In this question, we will focus on all crimes in the city of Chicago. Subsequently, we will look at specific crimes.

One way to look at how all crimes in Chicago have changed over time would be simply to plot number of crimes over time. However, this treats all crimes as being equally weighted, when in reality they are not. The RAND Corporation published a seminal study in which a team of economists established the costs of certain types of crime to society. We can create a crime index using the cost of crime to society, and examine how the cost of crime has changed over time.

The cost of crime numbers only exist for index crimes and are as follows:

First, we import the data. Since we are looking citywide, we will keep only the FBI Code and the Year columns at this point. In addition, we will only keep index crimes in this dataset. One thing I’ve noticed is that there do appear to be duplicates, so we keep only unique crimes by ID.

suppressPackageStartupMessages(library(readr))
suppressPackageStartupMessages(library(ggplot2))
## Warning: package 'ggplot2' was built under R version 3.3.3
suppressPackageStartupMessages(library(ggthemes))
## Warning: package 'ggthemes' was built under R version 3.3.3
suppressPackageStartupMessages(library(raster))
## Warning: package 'raster' was built under R version 3.3.3
## Warning: package 'sp' was built under R version 3.3.3
suppressPackageStartupMessages(library(dplyr))
## Warning: package 'dplyr' was built under R version 3.3.3
path = 'C:/Users/tdneumann/Downloads/ELUCD/Clean Data'
crime_files = list.files(path = path)
cost_of_crime = data.frame(FBI.Code = c('01A','01B','02','03','04A','04B','05','06','07'),
                           Cost_of_crime = c(rep(8649216, 2),217866,67266,rep(87238, 2), 13096, 2139, 9079))



cost_of_crime_summary = data.frame()

for (i in 1:length(crime_files)) {
  tmp = read_csv(paste(path, crime_files[i], sep = '/')) %>%
    distinct(ID,  .keep_all = TRUE) %>%
    select(`FBI Code`, Year) %>%
    rename(FBI.Code = `FBI Code`) %>%
    inner_join(cost_of_crime) %>%
    group_by(Year) %>%
    summarise(cost_of_crime_yearly = sum(Cost_of_crime)) %>%
    ungroup()
    
  cost_of_crime_summary = rbind(cost_of_crime_summary, tmp)
}
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   .default = col_character(),
##   X1 = col_integer(),
##   ID = col_integer(),
##   Beat = col_integer(),
##   District = col_double(),
##   Ward = col_double(),
##   `Community Area` = col_double(),
##   `X Coordinate` = col_double(),
##   `Y Coordinate` = col_double(),
##   Year = col_integer(),
##   Latitude = col_double(),
##   Longitude = col_double()
## )
## See spec(...) for full column specifications.
## Joining, by = "FBI.Code"
## Warning: Column `FBI.Code` joining character vector and factor, coercing
## into character vector
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   .default = col_character(),
##   X1 = col_integer(),
##   ID = col_integer(),
##   Beat = col_integer(),
##   District = col_double(),
##   Ward = col_double(),
##   `Community Area` = col_double(),
##   `X Coordinate` = col_double(),
##   `Y Coordinate` = col_double(),
##   Year = col_integer(),
##   Latitude = col_double(),
##   Longitude = col_double()
## )
## See spec(...) for full column specifications.
## Joining, by = "FBI.Code"
## Warning: Column `FBI.Code` joining character vector and factor, coercing
## into character vector

## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   .default = col_character(),
##   X1 = col_integer(),
##   ID = col_integer(),
##   Beat = col_integer(),
##   District = col_double(),
##   Ward = col_double(),
##   `Community Area` = col_double(),
##   `X Coordinate` = col_double(),
##   `Y Coordinate` = col_double(),
##   Year = col_integer(),
##   Latitude = col_double(),
##   Longitude = col_double()
## )
## See spec(...) for full column specifications.
## Joining, by = "FBI.Code"
## Warning: Column `FBI.Code` joining character vector and factor, coercing
## into character vector

## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   .default = col_character(),
##   X1 = col_integer(),
##   ID = col_integer(),
##   Beat = col_integer(),
##   District = col_double(),
##   Ward = col_double(),
##   `Community Area` = col_double(),
##   `X Coordinate` = col_double(),
##   `Y Coordinate` = col_double(),
##   Year = col_integer(),
##   Latitude = col_double(),
##   Longitude = col_double()
## )
## See spec(...) for full column specifications.
## Joining, by = "FBI.Code"
## Warning: Column `FBI.Code` joining character vector and factor, coercing
## into character vector
g = ggplot(subset(cost_of_crime_summary, Year < 2017)) +
  geom_line(aes(x = Year, y = cost_of_crime_yearly), size = 1.7, colour = 'indianred4') +
  scale_x_continuous(breaks = 2001:2016) +
  scale_y_continuous(breaks = seq(5, 12, by = 1)*10^9, labels = scales::comma) +
  ylab('Cost of Crime Anually (in $)') +
  ggtitle('Cost of Index Crimes to Society from 2001 - 2016\nin the City of Chicago') +
  theme_gdocs()

g

This chart shows interesting information regarding how the cost of crime has changed over time in the City of Chicago. We only have a single partial month of data for 2017, so we can ignore the low value for 2017.

We can see interesting patterns between these cost of crime numbers and leadership changes at CPD. 2001 and 2002 represent two of the worst years in the series, and unsurprisingly, at the start of 2003, a leadership change was made from Supt Hillard to Supt Cline. Supt Cline presided over much lower years than Hillard, but saw a spike in 2006. He resigned at the beginning of 2007, and Supt Weis began his tenure the following year. Weis’ crime fighting strategies generally continued the downward trend of cost of crime, and he resigned in early 2011. Garry McCarthy from New York stepped into the role in mid 2011, and presided over the three lowest cost of crime years since 2001. However, despite his strong record reducing crime, he was asked to resign by Mayor Emmanuel in Dec 2015 due to the cover-up of the fatal shooting of an unarmed citizen, Laquan McDonald. 2016 then saw the new Supt Eddie Johnson take office, along with a significant uptick in the cost of crime, due largely to an increase in homicides and shootings.

Now, let’s examine 120 x 120 grid – called a raster layer – and aggregate all of the crime yearly into these squares. We can then create a heatmap from these squares to show where the cost of crime is the most, and more interestingly, how it has changed, over time.

suppressMessages(library(rgdal))
suppressMessages(library(sp))

#### Shape file for City of Chicago Police districts
districts = readOGR('C:/Users/tdneumann/Documents/Districts')
## OGR data source with driver: ESRI Shapefile 
## Source: "C:/Users/tdneumann/Documents/Districts", layer: "Districts"
## with 25 features
## It has 2 fields
cost_of_crime_location = data.frame()


### Create file for lat / long and cost of crime since 2001
for (i in 1:length(crime_files)) {
  
  tmp = read_csv(paste(path, crime_files[i], sep = '/')) %>%
    distinct(ID,  .keep_all = TRUE) %>%
    select(`FBI Code`, Year, Latitude, Longitude) %>%
    rename(FBI.Code = `FBI Code`) %>%
    inner_join(cost_of_crime)
    
  cost_of_crime_location = rbind(cost_of_crime_location, tmp)
}
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   .default = col_character(),
##   X1 = col_integer(),
##   ID = col_integer(),
##   Beat = col_integer(),
##   District = col_double(),
##   Ward = col_double(),
##   `Community Area` = col_double(),
##   `X Coordinate` = col_double(),
##   `Y Coordinate` = col_double(),
##   Year = col_integer(),
##   Latitude = col_double(),
##   Longitude = col_double()
## )
## See spec(...) for full column specifications.
## Joining, by = "FBI.Code"
## Warning: Column `FBI.Code` joining character vector and factor, coercing
## into character vector

## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   .default = col_character(),
##   X1 = col_integer(),
##   ID = col_integer(),
##   Beat = col_integer(),
##   District = col_double(),
##   Ward = col_double(),
##   `Community Area` = col_double(),
##   `X Coordinate` = col_double(),
##   `Y Coordinate` = col_double(),
##   Year = col_integer(),
##   Latitude = col_double(),
##   Longitude = col_double()
## )
## See spec(...) for full column specifications.
## Joining, by = "FBI.Code"
## Warning: Column `FBI.Code` joining character vector and factor, coercing
## into character vector

## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   .default = col_character(),
##   X1 = col_integer(),
##   ID = col_integer(),
##   Beat = col_integer(),
##   District = col_double(),
##   Ward = col_double(),
##   `Community Area` = col_double(),
##   `X Coordinate` = col_double(),
##   `Y Coordinate` = col_double(),
##   Year = col_integer(),
##   Latitude = col_double(),
##   Longitude = col_double()
## )
## See spec(...) for full column specifications.
## Joining, by = "FBI.Code"
## Warning: Column `FBI.Code` joining character vector and factor, coercing
## into character vector

## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   .default = col_character(),
##   X1 = col_integer(),
##   ID = col_integer(),
##   Beat = col_integer(),
##   District = col_double(),
##   Ward = col_double(),
##   `Community Area` = col_double(),
##   `X Coordinate` = col_double(),
##   `Y Coordinate` = col_double(),
##   Year = col_integer(),
##   Latitude = col_double(),
##   Longitude = col_double()
## )
## See spec(...) for full column specifications.
## Joining, by = "FBI.Code"
## Warning: Column `FBI.Code` joining character vector and factor, coercing
## into character vector
# Remove points missing lat / long
cost_of_crime_location = subset(cost_of_crime_location, !is.na(Latitude))

years = 2001:2016

for (i in 1:length(years)) {
  
  cost_of_crime_location_sub = subset(cost_of_crime_location, Year == years[i])
  cost_of_crime_sp = SpatialPointsDataFrame(coords = cost_of_crime_location_sub[,c('Longitude', 'Latitude')], 
                                            data = cost_of_crime_location_sub, 
                                            proj4string = districts@proj4string)
    r = raster(ncol=120, nrow=120)
    extent(r) = extent(districts)
    proj4string(r) = proj4string(districts)
    r2 = rasterize(cost_of_crime_sp, r, cost_of_crime_sp$Cost_of_crime, fun=sum)
    r2[r2 >= 1.8*10^7] = 1.8*10^7
    plot(r2, main = years[i], col = colorRampPalette(c("lightgrey", "orange", "darkred"))(255))
    
}

Above, we see that the cost of crime has changed spatially quite significantly since 2001. I capped the value at 15 million dollars so that the scale would remain consistent from year to year. In 2001, it appears to be distributed more evenly throughout the city. As we saw with the previous figure, this was also the most costly year on record for crime. From the map, we can see significant portions of roughly 2x2 block areas of the city experiencing at least $15,000,000 in societal costs of crime. We also see that the near north side to the far north side had high levels of crime until approximately 2010, when this seems to dissipate. This could be due in part to the demolition of Cabrini Green, a housing project on the near north side notorious for spawning criminal activity and organized crime.

As the years progress, especially after 2010, we see generally see the crime move further west and south. This is certainly the case from 2014 to 2016, and in 2016, we see the highest cost of crime coming from a small concentrated area on the west side of the city. This area saw a dramatic uptick in their homicides from the previous year.

If I had more time, I would like to make the previous plots into an animation so that we could clearly see how the location of crime changes over time. I would also work on the color palette so that changes in the scale could be seen more clearly.

In part 2, we will examine how specific crimes have changed over time.

How have certain crimes changed over time?

In this section, we will take a brief look at how certain types of crimes have changed over time. Again, we will only focus on index crimes.

crime_type_summary = data.frame()


### Create file for lat / long and cost of crime since 2001
for (i in 1:length(crime_files)) {
  
  tmp = read_csv(paste(path, crime_files[i], sep = '/')) %>%
    distinct(ID,  .keep_all = TRUE) %>%
    select(`FBI Code`, Year) %>%
    rename(FBI.Code = `FBI Code`) %>%
    filter(FBI.Code %in% c('01A','01B','02','03','04A','04B','05','06','07')) %>%
    group_by(Year, FBI.Code) %>% 
    summarise(count_crime = n()) %>%
    ungroup()

  crime_type_summary = rbind(crime_type_summary, tmp)
}
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   .default = col_character(),
##   X1 = col_integer(),
##   ID = col_integer(),
##   Beat = col_integer(),
##   District = col_double(),
##   Ward = col_double(),
##   `Community Area` = col_double(),
##   `X Coordinate` = col_double(),
##   `Y Coordinate` = col_double(),
##   Year = col_integer(),
##   Latitude = col_double(),
##   Longitude = col_double()
## )
## See spec(...) for full column specifications.
## Warning: package 'bindrcpp' was built under R version 3.3.3

## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   .default = col_character(),
##   X1 = col_integer(),
##   ID = col_integer(),
##   Beat = col_integer(),
##   District = col_double(),
##   Ward = col_double(),
##   `Community Area` = col_double(),
##   `X Coordinate` = col_double(),
##   `Y Coordinate` = col_double(),
##   Year = col_integer(),
##   Latitude = col_double(),
##   Longitude = col_double()
## )
## See spec(...) for full column specifications.
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   .default = col_character(),
##   X1 = col_integer(),
##   ID = col_integer(),
##   Beat = col_integer(),
##   District = col_double(),
##   Ward = col_double(),
##   `Community Area` = col_double(),
##   `X Coordinate` = col_double(),
##   `Y Coordinate` = col_double(),
##   Year = col_integer(),
##   Latitude = col_double(),
##   Longitude = col_double()
## )
## See spec(...) for full column specifications.
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   .default = col_character(),
##   X1 = col_integer(),
##   ID = col_integer(),
##   Beat = col_integer(),
##   District = col_double(),
##   Ward = col_double(),
##   `Community Area` = col_double(),
##   `X Coordinate` = col_double(),
##   `Y Coordinate` = col_double(),
##   Year = col_integer(),
##   Latitude = col_double(),
##   Longitude = col_double()
## )
## See spec(...) for full column specifications.
#remove 01B, there are too few every year

crime_type_summary = subset(crime_type_summary, FBI.Code != '01B' & Year < 2017)

ggplot(crime_type_summary) +
  geom_line(aes(x = Year, y = count_crime, colour = FBI.Code), size = 1.7) +
  facet_wrap(~FBI.Code, nrow = 4,scales = 'free_y') +
  theme_gdocs() + 
  xlab('Year') +
  ylab('Number of Crimes per Year') +
  ggtitle('Yearly Index Crimes in the City of Chicago\n2001 - 2016')

Here we see a very interesting pattern in the data. For all but one category – homicides – the data seems to show a consistent pattern. We see a generally decreasing crime pattern from 2001 - 2004, a slight increase in 2005, and then a significant spike in 2006, followed by a small decrease in 2007, and a peak in 2008. Levels remain quite high from 2008 to 2011, until Supt McCarthy takes office, and levels remain consistently low.

Homicides are a different story. Homicides are catastrophically high in 2001 and 2002, but after this period they decrease significantly and remain between 400 and 600 per year until 2016, when they hit a high of 769, just 10 fewer than in 2002, reversing a 15 year trend.

Can you forecast how crime is likely to change in the future?

I have spent significant time working on forecasting crime into the future for the Chicago Police Department. In fact, I built a software tool that produces monthly forecasts for every crime in every district for the city. I use prophet - an automatic forecasting package developed by Facebook’s internal R&D team – and a cross-validation method developed by Hyndman to create models that can forecast crime accurately and produce confidence intervals. Below, I highlight the cross validation framework for choosing the best parameters. We will use Hyndman’s Evaluation on a rolling forecasting origin method, illustrated by the graph below. The blue points represent training periods, and the red point represents validation periods. We roll over several validation periods, adding a new training period each time, and averaging our evaluation metric across all of the validation periods. Hyndman’s cv framework

Let’s create a forecast for monthly robberies in the city of Chicago using this method. We will use 12 cv periods - each month in 2014, then forecast into 2015 to seen our forecast error.

suppressPackageStartupMessages(library(zoo))
suppressPackageStartupMessages(library(prophet))
## Warning: package 'Rcpp' was built under R version 3.3.3
prophet_frame = data.frame()

for (i in 1:length(crime_files)) {
  
  tmp = read_csv(paste(path, crime_files[i], sep = '/')) %>%
    distinct(ID,  .keep_all = TRUE) %>%
    select(`FBI Code`, Date) %>%
    rename(FBI.Code = `FBI Code`) %>%
    filter(FBI.Code == '03') %>%
    mutate(month_year = as.Date(as.yearmon(as.POSIXct(Date, format = '%m/%d/%Y %I:%M:%S %p')))) %>%
    group_by(month_year, FBI.Code) %>% 
    summarise(count_crime = n()) %>%
    ungroup() %>%
    select(month_year, count_crime) %>%
    rename(ds = month_year, y = count_crime)

  prophet_frame = rbind(prophet_frame, tmp)
}
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   .default = col_character(),
##   X1 = col_integer(),
##   ID = col_integer(),
##   Beat = col_integer(),
##   District = col_double(),
##   Ward = col_double(),
##   `Community Area` = col_double(),
##   `X Coordinate` = col_double(),
##   `Y Coordinate` = col_double(),
##   Year = col_integer(),
##   Latitude = col_double(),
##   Longitude = col_double()
## )
## See spec(...) for full column specifications.
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   .default = col_character(),
##   X1 = col_integer(),
##   ID = col_integer(),
##   Beat = col_integer(),
##   District = col_double(),
##   Ward = col_double(),
##   `Community Area` = col_double(),
##   `X Coordinate` = col_double(),
##   `Y Coordinate` = col_double(),
##   Year = col_integer(),
##   Latitude = col_double(),
##   Longitude = col_double()
## )
## See spec(...) for full column specifications.
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   .default = col_character(),
##   X1 = col_integer(),
##   ID = col_integer(),
##   Beat = col_integer(),
##   District = col_double(),
##   Ward = col_double(),
##   `Community Area` = col_double(),
##   `X Coordinate` = col_double(),
##   `Y Coordinate` = col_double(),
##   Year = col_integer(),
##   Latitude = col_double(),
##   Longitude = col_double()
## )
## See spec(...) for full column specifications.
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   .default = col_character(),
##   X1 = col_integer(),
##   ID = col_integer(),
##   Beat = col_integer(),
##   District = col_double(),
##   Ward = col_double(),
##   `Community Area` = col_double(),
##   `X Coordinate` = col_double(),
##   `Y Coordinate` = col_double(),
##   Year = col_integer(),
##   Latitude = col_double(),
##   Longitude = col_double()
## )
## See spec(...) for full column specifications.
set.seed(12345)
origin_month  = as.Date('2015-07-01' )
end_month     = as.Date('2016-07-01') 
date_seq      = seq(origin_month, end_month, by = 'month')

cv_set  = prophet_frame %>%
          filter(ds <= end_month)

# 21 random points -- the first two are the base parameters

rand_search_grid =  data.frame( 
  changepoint_prior_scale = c(0.05, sort(runif(20, 0.01, 0.1))),
  seasonality_prior_scale = c(10, sort(sample(c(runif(10, 0.01, 0.05), runif(10, 1, 10)), 10, replace = F)),
                              sort(sample(c(runif(10, 0.01, 0.05), runif(10, 1, 10)), 10, replace = F))),
  Value = rep(0, 21)
)


  # Search best parameters
for (i in 1:nrow(rand_search_grid)) {
  parameters = rand_search_grid[i, ]
  error = c()
  for (d in date_seq) {
    train = subset(cv_set, ds < as.Date(d, origin = '1970-01-01'))
    test  = subset(cv_set, ds == as.Date(d, origin = '1970-01-01'))
    
    m = prophet(train, growth = 'linear',
                 seasonality.prior.scale = parameters$seasonality_prior_scale, 
                 changepoint.prior.scale = parameters$changepoint_prior_scale,
                 weekly.seasonality = F,
                 daily.seasonality = F)
    
    future = make_future_dataframe(m, periods = 1, freq = 'month')
    forecast = predict(m, future)
    forecast$ds = as.Date(forecast$ds)
    
    error_d = forecast::accuracy(forecast[forecast$ds %in% test$ds, 'yhat'], test$y)[ , 'MAPE']
    error = c(error, error_d)
  }
  mean_error = mean(error)
  print(mean_error)
  rand_search_grid$Value[i] = -mean_error
}
## Initial log joint probability = -4.45569
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.14737
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.13858
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.12761
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.12
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.2917
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.17486
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.16596
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.60358
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.37843
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.38592
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.13899
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.14262
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## [1] 16.52868
## Initial log joint probability = -4.45569
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.14737
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.13858
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.12761
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.12
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.2917
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.17486
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.16596
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.60358
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.37843
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.38592
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.13899
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.14262
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## [1] 11.5348
## Initial log joint probability = -4.45569
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.14737
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.13858
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.12761
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.12
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.2917
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.17486
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.16596
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.60358
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.37843
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.38592
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.13899
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.14262
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## [1] 13.33764
## Initial log joint probability = -4.45569
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.14737
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.13858
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.12761
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.12
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.2917
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.17486
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.16596
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.60358
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.37843
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.38592
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.13899
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.14262
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## [1] 15.10911
## Initial log joint probability = -4.45569
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.14737
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.13858
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.12761
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.12
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.2917
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.17486
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.16596
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.60358
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.37843
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.38592
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.13899
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.14262
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## [1] 15.01939
## Initial log joint probability = -4.45569
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.14737
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.13858
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.12761
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.12
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.2917
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.17486
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.16596
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.60358
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.37843
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.38592
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.13899
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.14262
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## [1] 16.56946
## Initial log joint probability = -4.45569
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.14737
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.13858
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.12761
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.12
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.2917
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.17486
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.16596
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.60358
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.37843
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.38592
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.13899
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.14262
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## [1] 16.39542
## Initial log joint probability = -4.45569
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.14737
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.13858
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.12761
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.12
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.2917
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.17486
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.16596
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.60358
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.37843
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.38592
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.13899
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.14262
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## [1] 16.35464
## Initial log joint probability = -4.45569
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.14737
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.13858
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.12761
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.12
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.2917
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.17486
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.16596
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.60358
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.37843
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.38592
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.13899
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.14262
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## [1] 16.76207
## Initial log joint probability = -4.45569
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.14737
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.13858
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.12761
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.12
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.2917
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.17486
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.16596
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.60358
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.37843
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.38592
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.13899
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.14262
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## [1] 16.26162
## Initial log joint probability = -4.45569
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.14737
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.13858
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.12761
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.12
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.2917
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.17486
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.16596
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.60358
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.37843
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.38592
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.13899
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.14262
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## [1] 17.62574
## Initial log joint probability = -4.45569
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.14737
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.13858
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.12761
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.12
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.2917
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.17486
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.16596
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.60358
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.37843
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.38592
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.13899
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.14262
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## [1] 12.80742
## Initial log joint probability = -4.45569
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.14737
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.13858
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.12761
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.12
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.2917
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.17486
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.16596
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.60358
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.37843
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.38592
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.13899
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.14262
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## [1] 14.61579
## Initial log joint probability = -4.45569
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.14737
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.13858
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.12761
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.12
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.2917
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.17486
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.16596
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.60358
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.37843
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.38592
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.13899
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.14262
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## [1] 15.28808
## Initial log joint probability = -4.45569
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.14737
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.13858
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.12761
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.12
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.2917
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.17486
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.16596
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.60358
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.37843
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.38592
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.13899
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.14262
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## [1] 16.09053
## Initial log joint probability = -4.45569
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.14737
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.13858
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.12761
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.12
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.2917
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.17486
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.16596
## Optimization terminated normally: 
##   Convergence detected: relative change in objective function was below tolerance
## Initial log joint probability = -4.60358
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.37843
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.38592
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.13899
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.14262
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## [1] 16.33554
## Initial log joint probability = -4.45569
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.14737
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.13858
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.12761
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.12
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.2917
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.17486
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.16596
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.60358
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.37843
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.38592
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.13899
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.14262
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## [1] 16.39766
## Initial log joint probability = -4.45569
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.14737
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.13858
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.12761
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.12
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.2917
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.17486
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.16596
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.60358
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.37843
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.38592
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.13899
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.14262
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## [1] 19.35951
## Initial log joint probability = -4.45569
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.14737
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.13858
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.12761
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.12
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.2917
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.17486
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.16596
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.60358
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.37843
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.38592
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.13899
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.14262
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## [1] 19.563
## Initial log joint probability = -4.45569
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.14737
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.13858
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.12761
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.12
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.2917
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.17486
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.16596
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.60358
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.37843
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.38592
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.13899
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.14262
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## [1] 19.38804
## Initial log joint probability = -4.45569
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.14737
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.13858
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.12761
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.12
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.2917
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.17486
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.16596
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.60358
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.37843
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
## Initial log joint probability = -4.38592
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.13899
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## Initial log joint probability = -4.14262
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance
## [1] 19.75645
best_cv_value = -1*rand_search_grid[which.max(rand_search_grid$Value),'Value']
rand_search_grid = arrange(rand_search_grid, -Value)
head(rand_search_grid)
##   changepoint_prior_scale seasonality_prior_scale     Value
## 1              0.01010229              0.01905869 -11.53480
## 2              0.05162452              0.01946322 -12.80742
## 3              0.01310819              0.02814912 -13.33764
## 4              0.05583019              0.02742122 -14.61579
## 5              0.02497346              0.03794175 -15.01939
## 6              0.02371361              0.03578171 -15.10911
best_params = rand_search_grid[1,]

print('MAPE with base params:')
## [1] "MAPE with base params:"
print(subset(rand_search_grid, changepoint_prior_scale == 0.05)$Value*-1)
## [1] 16.52868
print('MAPE with tuned params:')
## [1] "MAPE with tuned params:"
print(best_params$Value*-1)
## [1] 11.5348

Now that we have the best parameters established, we can refit with all of the training data and test on our holdout set, all months from 2016.

test = subset(prophet_frame, ds >= end_month & ds < as.Date('2017-01-01'))

m = prophet(cv_set, growth = 'linear',
                 seasonality.prior.scale = best_params$seasonality_prior_scale, 
                 changepoint.prior.scale = best_params$changepoint_prior_scale)
## Disabling weekly seasonality. Run prophet with weekly.seasonality=TRUE to override this.
## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.
## Initial log joint probability = -4.4296
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
future = make_future_dataframe(m, periods = 5, freq = 'month')
    
forecast = predict(m, future)
forecast$ds = as.Date(forecast$ds)

p = ggplot() + 
    geom_point(data = cv_set, aes(x = as.Date(ds), y = y), size = 1) +
    geom_line(data = forecast, aes(x = as.Date(ds), y = yhat), color = "#0072B2", size = 1) +
    geom_ribbon(data = forecast, aes(x = as.Date(ds), ymin = yhat_lower, ymax = yhat_upper), fill = "#0072B2", alpha = 0.3) +
    geom_point(data = test, aes(x = as.Date(ds), y = y), size = 1, color = '#4daf4a') 
p

We see that our holdout set (2016) seems to be higher than predicted, but well within the 80% confidence interval month over month.

It would be easy to functionalize the above code so that it could be run on any crime for any geography. Additionally, since I implemented random grid search, it would make sense to keep the best parameters from the previous run in a database and add them to a new set of random parameters every time the model is retrained (likely monthly)

Here’s a use of these forecasts in practice.

In 2016, The 006th district – a largely middle-class African-American neighborhood on Chicago’s far south side – experienced a dramatic increase in gun violence. This trend continued well into 2017. I was working as an analyst in the district at the time, and produced forecast models showing that shootings were predicted to be significantly higher for the summer of 2017 than for 2016, and the commander should try to get additional overtime or extra resoures from headquarters to help prevent this dramatic increase.