R, Python 분석과 프로그래밍의 친구 (by R Friend)

  1. Apache MADlib 소개

  2. web site에서 abalone 데이터셋 가져와서 table 만들기

---------------------------------

-- Linear Regression in Parallel 

-- using Apache MADlib

---------------------------------

-- Dataset for example: abalone dataset from the UC Irvine Machine Learning Repository

-- URL: http://archive.ics.uci.edu/ml/machine-learning-databases/abalone/abalone.data

-- Create an external web table

DROP EXTERNAL TABLE IF EXISTS abalone_external;

CREATE EXTERNAL WEB TABLE abalone_external(

sex text 

, length float8

, diameter float8

, height float8

, whole_weight float8

, shucked_weight float8

, viscera_weight float8

, shell_weight float8

, rings integer -- target variable to predict

) LOCATION('http://archive.ics.uci.edu/ml/machine-learning-databases/abalone/abalone.data') 

FORMAT 'CSV' 

(null as '?');

-- Create a table of abalone

DROP TABLE IF EXISTS abalone;

CREATE TABLE abalone AS 

SELECT * FROM abalone_external

DISTRIBUTED BY (sex);

-- Viewing data distribution

SELECT gp_segment_id, COUNT(*) AS row_cnt

FROM abalone

GROUP BY gp_segment_id;

-- Check data

SELECT * FROM abalone LIMIT 5;

SELECT sex, COUNT(*) FROM abalone GROUP BY sex; 

  3. 훈련, 검증 데이터 분할 (split abalone dataset into training and test set)

-- Train, Test set split

DROP TABLE IF EXISTS out_train, out_test;

SELECT madlib.train_test_split(

'abalone',  -- Source table

'out',      -- Output table

        0.8,        -- train_proportion

        NULL,       -- Default = 1 - train_proportion = 0.5

        'sex',      -- Strata definition

        'rings, shucked_weight, diameter', -- Columns to output

        FALSE,      -- Sample with replacement

        TRUE);      -- Separate output tables

-- Check

SELECT * FROM out_train LIMIT 5;

SELECT sex, count(*) FROM out_train GROUP BY sex;

SELECT sex, count(*) FROM out_test GROUP BY sex;

  4. 성별 그룹별 선형회귀모형 적합 

     (training linear regression models by 'sex' groups')

-- Linear Regression using MADlib

-- Train a regression model. 

DROP TABLE IF EXISTS abalone_linregr, abalone_linregr_summary;

SELECT madlib.linregr_train(

    'out_train'         -- table containing training data

    , 'abalone_linregr' -- table in which to save results

    , 'rings'           -- column containing dependent variable

    , 'ARRAY[1, shucked_weight, diameter]' -- features included in the model

    , 'sex'); -- create multiple output models (one for each value of sex)

  5. 성별 그룹별 선형회귀모형 회귀계수 조회 

     (select coefficients per variables by 'sex' groups)

-- Examine the resulting models

SELECT * FROM abalone_linregr ORDER BY sex;

-- unnest format

SELECT sex

, unnest(ARRAY['intercept', 'rings', 'diameter']) as attribute

, unnest(coef) as coefficient 

, unnest(std_err) as standard_error

, unnest(t_stats) as t_stat

, unnest(p_values) as pvalue

FROM abalone_linregr

ORDER BY sex;

  6. 성별 그룹별 선형회귀모형을 이용하여 예측 (prediction by 'sex' groups)

-- compare predicted value with actual with grouping

DROP TABLE IF EXISTS abalone_pred;

CREATE TABLE abalone_pred AS (

SELECT a.sex, a.shucked_weight, a.diameter, a.rings, 

madlib.linregr_predict(m.coef

, ARRAY[1, shucked_weight, diameter]

) as predict_val

, rings - madlib.linregr_predict(m.coef

, ARRAY[1, shucked_weight, diameter]

) as residual

FROM out_test a, abalone_linregr m

WHERE a.sex = m.sex) DISTRIBUTED BY (sex);

SELECT * FROM abalone_pred WHERE sex = 'F' LIMIT 10;

-- Model Performance Evaluation: Mean Squared Error

DROP TABLE IF EXISTS abalone_mse;

SELECT madlib.mean_squared_error(

'abalone_pred'   -- table_in

, 'abalone_mse'  -- table_out

, 'predict_val'  -- prediction_col

, 'rings'        -- observed_col

, 'sex');        -- grouping_cols

SELECT * FROM abalone_mse;

-- Mean absolute error

madlib.mean_abs_error(table_in, table_out, prediction_col, observed_col, grouping_cols)

-- Mean absolute percentage error

madlib.mean_abs_perc_error(table_in, table_out, prediction_col, observed_col, grouping_cols)

-- Mean percentage error

madlib.mean_perc_error(table_in, table_out, prediction_col, observed_col, grouping_cols)

-- Mean squared error

madlib.mean_squared_error(table_in, table_out, prediction_col, observed_col, grouping_cols)

-- R-squared

madlib.r2_score(table_in, table_out, prediction_col, observed_col, grouping_cols)

-- Adjusted R-squared

madlib.adjusted_r2_score(table_in, table_out, prediction_col, observed_col, num_predictors, training_size, grouping_cols)

----------------------------------

-- PL/R on Greenplum

-- : 5 ways to return PL/R results

----------------------------------

-- create an example dataset table

create schema test;

drop table if exists test.src_tbl;

create table test.src_tbl (

grp varchar(10) not null 

, x1 integer

, x2 float8

) distributed by (grp);

insert into test.src_tbl values 

('a', 1, 0.13)

, ('a', 2, 0.34)

, ('a', 3, 0.31)

, ('a', 4, 0.49)

, ('a', 5, 0.51)

, ('b', 1, 0.10)

, ('b', 2, 0.26)

, ('b', 3, 0.30)

, ('b', 4, 0.62)

, ('b', 5, 0.59);

select * from test.src_tbl;

  (1) returns float8[] : array 형태로 결과 반환 (그룹 당 1행)

------------------------------------------------

-- (1) returns an array using 'returns float8[]'

------------------------------------------------

-- Define PL/R UDF

drop function if exists test.plr_log_trans(integer[], float8[]);

create or replace function test.plr_log_trans(x1 integer[], x2 float8[])

returns float8[] -- returns as an array

as 

$$

x_new <- 2 * x1 + log(x2)

x_new <- round(x_new, 3)

return(x_new)

$$ language 'plr';

-- Execute PL/R UDF

select 

a.grp

, a.x1_agg

, a.x2_agg

, test.plr_log_trans(a.x1_agg, a.x2_agg) as x_new_agg

from (

select 

grp

, array_agg(x1::int) as x1_agg

, array_agg(x2::float8) as x2_agg

from test.src_tbl

group by grp

) a;

-- Display by rows using unnest() function

select 

a.grp

, unnest(a.x1_agg) as x1

, unnest(a.x2_agg) as x2

, unnest(test.plr_log_trans(a.x1_agg, a.x2_agg)) as x_new

from (

select 

grp

, array_agg(x1::int) as x1_agg

, array_agg(x2::float8) as x2_agg

from test.src_tbl

group by grp

) a;

  (2) returns setof float8 : 행 단위로 결과 반환 (관측치 당 1행, 그룹 당 여러 행)

-----------------------------------------

-- (2) returns rows using 'returns setof'

-----------------------------------------

-- Define PL/R UDF

drop function if exists test.plr_log_trans_2(int[], float8[]);

create or replace function test.plr_log_trans_2(x1 int[], x2 float8[])

returns setof float8 -- returns by each row

as 

$$

x_new <- 2 * x1 + log(x2)

x_new <- round(x_new, 3)

return(x_new)

$$ language 'plr';

-- Execute PL/R UDF

select 

a.grp

, unnest(a.x1_agg) as x1

, unnest(a.x2_agg) as x2

, test.plr_log_trans_2(a.x1_agg, a.x2_agg) as x_new

from (

select 

grp

, array_agg(x1::int) as x1_agg

, array_agg(x2::float8) as x2_agg

from test.src_tbl

group by grp

) a;

  (3) reterns setof composite_type : 행 단위로 composite type 에 맞게 결과 반환 

                                                (관측치 당 1행, 그룹 당 여러 행)

-------------------------------------------------------------

-- (3) returns rows using 'returns setof' with a composite type

-------------------------------------------------------------

-- Define composite data type

drop type test.plr_log_trans_type cascade;

create type test.plr_log_trans_type as (

x1 integer

, x2 float8

, x_new float8

);

-- Define PL/R UDF

drop function if exists test.plr_log_trans_3(int[], float8[]);

create or replace function test.plr_log_trans_3(x1 int[], x2 float8[])

returns setof test.plr_log_trans_type

as 

$$

x_new <- 2 * x1 + log(x2)

x_new <- round(x_new, 3)

result_df <- data.frame(x1, x2, x_new)

return(result_df)

$$ language 'plr';

-- Execute PL/R UDF

select 

a.grp

, (test.plr_log_trans_3(a.x1_agg, a.x2_agg)).*

from (

select 

grp

, array_agg(x1::int) as x1_agg

, array_agg(x2::float8) as x2_agg

from test.src_tbl

group by grp

) a;

  (4) returns table : 행 단위로 결과 반환 (관측치 당 1행, 그룹 당 여러 행)

-------------------------------------------

-- (4) returns rows using 'returns table'

-------------------------------------------

-- Define PL/R UDF

drop function if exists test.plr_log_trans_4(integer[], float8[]);

create or replace function test.plr_log_trans_4(x1 integer[], x2 float8[])

returns table (

x1 integer

, x2 float8

, x_new float8

)

as 

$$

x_new <- 2 * x1 + log(x2)

x_new <- round(x_new, 3)

result_df <- data.frame(x1, x2, x_new)

return(result_df)

$$ language 'plr';

-- Execute PL/R UDF

select 

a.grp

, (test.plr_log_trans_4(a.x1_agg, a.x2_agg)).*

from (

select 

grp

, array_agg(x1::int) as x1_agg

, array_agg(x2::float8) as x2_agg

from test.src_tbl

group by grp

) a;

  (5) returns bytea : byte array 로 결과 반환 (그룹 당 1행)

       --> unserialize 하는 PL/R 함수 추가로 필요

----------------------------------------------------------

-- (5) returns serialized byte array using 'returns bytea'

----------------------------------------------------------

-- Define PL/R UDF

drop function if exists test.plr_log_trans_5(integer[], float8[]);

create or replace function test.plr_log_trans_5(x1 integer[], x2 float8[])

returns bytea -- serialized as a byte array

as 

$$

x_new <- 2 * x1 + log(x2)

x_new <- round(x_new, 3)

serialized_df <- serialize(data.frame(x1, x2, x_new), NULL)

return(serialized_df)

$$ language 'plr';

-- Execute PL/R UDF

drop table if exists test.tbl_plr_log_trans_5;

create table test.tbl_plr_log_trans_5 as (

select 

a.grp

, test.plr_log_trans_5(a.x1_agg, a.x2_agg) as serialized_df

from (

select 

grp

, array_agg(x1::int) as x1_agg

, array_agg(x2::float8) as x2_agg

from test.src_tbl

group by grp

) a

) distributed by (grp);

select * from test.tbl_plr_log_trans_5;

-- Unserialize

-- Define composite data type

drop type test.plr_log_trans_type cascade;

create type test.plr_log_trans_type as (

x1 integer

, x2 float8

, x_new float8

);

-- Define PL/R UDF for reading a serialized PL/R results

drop function if exists test.plr_unserialize(bytea);

create or replace function test.plr_unserialize(

serialized_obj bytea

) returns setof test.plr_log_trans_type

as 

$$

unserialized_df <- data.frame(unserialize(serialized_obj))

return (unserialized_df)

$$ language 'plr';

-- Execute 'plr_unserialize' UDF

select 

grp

, (test.plr_unserialize(a.serialized_df)).*

from test.tbl_plr_log_trans_5 a;

--create schema and sample table

create schema test;

drop table if exists test.tbl;

create table test.tbl (

grp text not null

, x1 int

, x2 int

);

insert into test.tbl (grp, x1, x2) values 

('a', 1, 2)

, ('a', 2, 5)

, ('a', 3, 4)

, ('b', 1, 8)

, ('b', 2, 7)

, ('b', 3, 3);

select * from test.tbl;

  1. 열(column)을 기준으로 여러개의 행(row)을 그룹별로 array aggregation

-----------------------------------------------

-- (1) data preparation : array_agg() by column

-----------------------------------------------

drop table if exists test.array_agg_by_col;

create table test.array_agg_by_col as (

select 

grp

, array_agg(x1) as x1_agg

, array_agg(x2) as x2_agg

from test.tbl 

group by grp

) distributed by (grp);

select * from test.array_agg_by_col order by grp;

-- define PL/R UDF

drop function if exists test.plr_cor(int[], int[]);

create or replace function test.plr_cor(x1 int[], x2 int[]) 

returns float8

as

$$ 

corr_coef <- cor(x1, x2)

return (corr_coef)

$$ language 'plr';

-- execute PL/R UDF

select 

grp

, test.plr_cor(x1_agg, x2_agg) as corr_coef

from test.array_agg_by_col;

  2. 행과 열을 기준으로 그룹별로 2D (2-dimensional) array aggregation  

------------------------------------------------------

-- (2) data preparation : 2D array MADlib matrix_agg()

--     : only with same data types

------------------------------------------------------

drop table if exists test.tbl_matrix_agg;

create table test.tbl_matrix_agg as (

select 

grp

, madlib.matrix_agg(array[x1, x2]) as mat_agg

from test.tbl 

group by grp

) distributed by (grp);

select * from test.tbl_matrix_agg;

-- define PL/R UDF

drop function if exists test.plr_cor_2(float8[]);

create or replace function test.plr_cor_2(mat_agg float8[]) 

returns float8

as

$$

df <- data.frame(mat_agg)

colnames(df) <- c("x1", "x2")

corr_coef <- with(df, cor(x1, x2))

return (corr_coef)

$$ language 'plr';

-- execute PL/R UDF

select 

grp

, test.plr_cor_2(mat_agg) as corr_coef

from test.tbl_matrix_agg;

  3. 문자열로 string aggregation 한 후 PL/R 코드 안에서 R 로 데이터셋 변환하기

--------------------------------------------------

-- (3-a) data preparation : string_agg() by column

--------------------------------------------------

drop table if exists test.string_agg_by_col;

create table test.string_agg_by_col as (

select 

grp

, string_agg(x1::text, '|') as x1_str_agg

, string_agg(x2::text, '|') as x2_str_agg

from test.tbl 

group by grp

) distributed by (grp);

select * from test.string_agg_by_col order by grp;

-- define PL/R UDF : (3-a) using DataFrame

drop function if exists test.plr_cor_3(text, text);

create or replace function test.plr_cor_3(x1 text, x2 text) 

returns float8

as

$$ 

# make a temp DataFrame

df_tmp <- data.frame(x1, x2)

# split by delimiter and reshape it in a long format

split_func <- function(x){

options(stringsAsFactors = FALSE) # Not to read strings as factors

df_split <- as.data.frame(

do.call('cbind'

, strsplit(as.character(x)

                                     # set delimiter with yours

    , split="|"

    , fixed=T)))

return (df_split)

}

df <- data.frame(lapply(df_tmp, split_func))

colnames(df) <- c("x1", "x2") # set column names

# convert a data type from text to numeric

df <- data.frame(sapply(df, as.numeric))

# calculate correlation coefficients

corr_coef <- with(df, cor(x1, x2))

return (corr_coef)

$$ language 'plr';

-- execute PL/R UDF

select 

grp

, test.plr_cor_3(x1_str_agg, x2_str_agg) as corr_coef

from test.string_agg_by_col

order by grp asc;

--------------------------------------------------

-- (3-b) data preparation : string_agg() by column

--------------------------------------------------

drop table if exists test.string_agg_by_col;

create table test.string_agg_by_col as (

select 

grp

, string_agg(x1::text, '|') as x1_str_agg

, string_agg(x2::text, '|') as x2_str_agg

from test.tbl 

group by grp

) distributed by (grp);

select * from test.string_agg_by_col order by grp;

-- define PL/R UDF : (3-b) using data.table tstrsplit() function

drop function if exists test.plr_cor_4(text, text);

create or replace function test.plr_cor_4(x1 text, x2 text) 

returns float8

as

$$ 

library(data.table)

# make a temp DataTable

dt_tmp <- data.table(x1, x2)

# split by delimiter and reshape it in a long format

dt_split_func <- function(x){

dt_split <- data.table(tstrsplit(x, split="|", fixed=T))

return(dt_split)

}

df <- data.frame(lapply(dt_tmp, dt_split_func))

colnames(df) <- c("x1", "x2") # set column names

# convert a data type from text to numeric

df <- data.frame(sapply(df, as.numeric))

# calculate correlation coefficients

corr_coef <- with(df, cor(x1, x2))

return (corr_coef)

$$ language 'plr';

-- execute PL/R UDF

select 

grp

, test.plr_cor_4(x1_str_agg, x2_str_agg) as corr_coef

from test.string_agg_by_col;

  (1) abalone 데이터셋으로 테이블 만들기

---------------------------------

-- Linear Regression in Parallel 

-- using PL/R

---------------------------------

-- Dataset for example: abalone dataset from the UC Irvine Machine Learning Repository

-- url: http://archive.ics.uci.edu/ml/machine-learning-databases/abalone/abalone.data

-- Create an external web table

DROP EXTERNAL TABLE IF EXISTS abalone_external;

CREATE EXTERNAL WEB TABLE abalone_external(

sex text 

, length float8

, diameter float8

, height float8

, whole_weight float8

, shucked_weight float8

, viscera_weight float8

, shell_weight float8

, rings integer -- target variable to predict

) LOCATION('http://archive.ics.uci.edu/ml/machine-learning-databases/abalone/abalone.data') 

FORMAT 'CSV' 

(null as '?');

-- Create a table of abalone

DROP TABLE IF EXISTS abalone;

CREATE TABLE abalone AS 

SELECT * FROM abalone_external

DISTRIBUTED BY (sex);

-- Viewing data distribution

SELECT gp_segment_id, COUNT(*) AS row_cnt

FROM abalone

GROUP BY gp_segment_id;

-- Check data

SELECT * FROM abalone LIMIT 5;

  (2) Train, Test set 테이블 분할 (train, test set split)

---------------------------

-- Train, Test set split

---------------------------

-- Check the number of observations per sex group(F, I, M)
SELECT sex, COUNT(*) FROM abalone GROUP BY sex;

-- Train, Test set split

DROP TABLE IF EXISTS out_train, out_test;

SELECT madlib.train_test_split(

'abalone',    -- Source table

'out',     -- Output table

    0.8,       -- train_proportion

    NULL,      -- Default = 1 - train_proportion = 0.5

    'sex', -- Strata definition

    'rings, shucked_weight, diameter', -- Columns to output

    FALSE,      -- Sample with replacement

    TRUE);     -- Separate output tables

SELECT * FROM out_train LIMIT 5;

SELECT sex, count(*) FROM out_train GROUP BY sex;

SELECT sex, count(*) FROM out_test GROUP BY sex;

  (3) array aggregation 하여 PL/R에서 사용할 데이터셋 준비하기

-- Data Preparation

-- : array aggregation using array_agg()

DROP TABLE IF EXISTS abalone_array;

CREATE TABLE abalone_array AS 

SELECT

sex::text -- group

, array_agg(rings::float8) as rings             -- y

, array_agg(shucked_weight::float8) as s_weight -- x1

, array_agg(diameter::float8) as diameter       -- x2

FROM out_train

GROUP BY sex

DISTRIBUTED BY (sex);

SELECT * FROM abalone_array;

  (4) 선형회귀모형 적합하는 PL/R 사용자 정의 함수 정의하기 (Define PL/R UDF)

----------------------------------------------------------------

-- (4-1) PL/R : Linear Regression Model's Coefficients --> Predict

----------------------------------------------------------------

-- Return Types

DROP TYPE IF EXISTS lm_abalone_type CASCADE;

CREATE TYPE lm_abalone_type AS (

variable text

, coef_est float

, std_error float

, t_stat float

, p_value float

);

-- PL/R User Defined Function

DROP FUNCTION IF EXISTS plr_lm_train(float8[], float8[], float8[]);

CREATE OR REPLACE FUNCTION plr_lm_train(

rings float8[]

, s_weight float8[]

, diameter float8[]

) RETURNS SETOF lm_abalone_type AS

$$

m1 <- lm(rings ~ s_weight + diameter)

m1_s <- summary(m1)$coef

temp_m1 <- data.frame(rownames(m1_s), m1_s)

return(temp_m1)

$$

LANGUAGE 'plr';

 (5) PL/R 실행하기 (execute PL/R UDF in Greenplum DB)

-- Execution of Linear Regression PL/R 

DROP TABLE IF EXISTS lm_abalone_model_coef;

CREATE TABLE lm_abalone_model_coef AS (

SELECT sex, (plr_lm_train(rings, s_weight, diameter)).* 

FROM abalone_array

) DISTRIBUTED BY (sex);

SELECT * FROM lm_abalone_model_coef;

  (6) 적합한 선형회귀모형을 사용해 test set에 대해 예측하기 (prediction on test set)

------------------------------------------------

-- Prediction and Model Evaluation for Test set

------------------------------------------------

-- Preparation of test set in aggregated array

DROP TABLE IF EXISTS test_array;

CREATE TABLE test_array AS 

SELECT

sex::text

, array_agg(rings::float8) as rings             -- y

, array_agg(shucked_weight::float8) as s_weight -- x1

, array_agg(diameter::float8) as diameter       --x2

FROM out_test

GROUP BY sex

DISTRIBUTED BY (sex);

SELECT * FROM test_array;

-- Define composite data type for predicted return values

DROP TYPE IF EXISTS lm_predicted_type CASCADE;

CREATE TYPE lm_predicted_type AS (

id int

, actual float

, pred float

);

-- Define PL/R UDF of prediction per groups using linear regression coefficients

DROP FUNCTION IF EXISTS plr_lm_coef_predict(float8[], float8[], float8[], float8[]);

CREATE OR REPLACE FUNCTION plr_lm_coef_predict(

rings float8[]

, s_weight float8[]

, diameter float8[]

, coef_est float8[]

) RETURNS SETOF lm_predicted_type AS

$$

actual <- rings # y

intercept <- 1

X <- cbind(intercept, s_weight, diameter) # X matrix

coef_est <- matrix(coef_est) # coefficients matrix

predicted <- X %*% coef_est  # matrix multiplication

df_actual_pred <- data.frame(actual, predicted)

id <- as.numeric(rownames(df_actual_pred))

return(data.frame(id, df_actual_pred))

$$

LANGUAGE 'plr';

-- Execute PL/R Prediction UDF

DROP TABLE IF EXISTS out_coef_predict;

CREATE TABLE out_coef_predict AS (

SELECT sex, (plr_lm_coef_predict(c.rings, c.s_weight, c.diameter, c.coef_est)).*

FROM (

SELECT a.*, b.coef_est

FROM test_array a, 

(SELECT sex, array_agg(coef_est) AS coef_est FROM lm_abalone_model_coef GROUP BY sex) b

WHERE a.sex = b.sex

) c

) DISTRIBUTED BY (sex);

-- Compare 'actual' vs. 'predicted'

SELECT * FROM out_coef_predict WHERE sex = 'F' ORDER BY sex, id LIMIT 10;

  (7) 회귀모형 자체를 Serialize 해서 DB에 저장하고, Unserialize 해서 예측하기

--------------------------------------------------------------------

-- (2) PL/R : Linear Model --> Serialize --> Deserialize --> Predict

--------------------------------------------------------------------

-- PL/R User Defined Function Definition

DROP FUNCTION IF EXISTS plr_lm_model(float8[], float8[], float8[]);

CREATE OR REPLACE FUNCTION plr_lm_model(

    rings float8[]

, s_weight float8[]

, diameter float8[]

) RETURNS bytea -- serialized model as a byte array

AS

$$

lr_model <- lm(rings ~ s_weight + diameter)

return (serialize(lr_model, NULL))

$$

LANGUAGE 'plr';

-- Execution of Linear Regression PL/R 

DROP TABLE IF EXISTS lm_abalone_model;

CREATE TABLE lm_abalone_model AS (

SELECT sex, plr_lm_model(rings, s_weight, diameter) AS serialized_model

FROM abalone_array

) DISTRIBUTED BY (sex);